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This  report  summarizes  the  research  activities  at  the  Information  Systems 
Laboratory,  Stanford  University,  for  the  Distributed  Sensor  Network  project  dur¬ 
ing  the  past  two  years.  The  objective  of  this  research  effort  was  to  develop  new 
signal  processing  techniques  for  locating,  tracking  and  classifying  multiple 
sources  and  receivers  in  a  DSN  system.  The  major  accomplishments  presented 
are  in  the  areas  of:  source  location  techniques,  estimation  algorithms  for  DSN, 
and  the  development  of  a  signal  processing  workbench.  A  comprehensive  over¬ 
view  of  our  research  results  is  presented  with  emphasis  on  their  applications  to 
various  facets  of  the  DSN  problem. 
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1.  INTRODUCTION 

This  report  summarizes  the  research  activities  in  the  Information  Systems 
Laboratory  (ISL)  at  Stanford  University  for  the  Distributed  Sensor  Network 
(DSN)  project.  Our  effort  is  part  of  a  DARPA  sponsored  program  with  several 
participant  contractors.  This  report  covers  the  period  16  December  1978  to  15 
December  1980. 

Our  research  effort  during  this  reporting  period  concentrated  on  three 
principal  areas: 

•  source  location  techniques 

•  estimation  algorithms  for  DSN 

•  development  of  a  signal  processing  workbench. 

This  work  continues  and  enhances  the  research  performed  during  the  first  year 
of  the  project,  as  reported  in  [DSN].  Our  major  accomplishments  are  briefly 
summarized  below. 

Sburce  location  Techniques 

Our  original  formulation  of  the  problem  of  estimating  source  location  from 
time  difference-of-arrival,  as  a  (conditionally)  linear  problem,  proved  to  be  very 
useful.  It  enabled  us  to  develop  for  the  first  time  optimal  location  estimation 
algorithms  and  to  perform  a  detailed  analysis  of  the  identifiability  issue.  The 
results  of  this  analysis  were  new  insights  into  the  sensor  placement  problem 
(Le.,  what  is  the  best  configuration  for  a  given  set  of  sensors  from  the  standpoint 
of  location  estimation  accuracy)  and  the  development  of  efficient  location  esti¬ 
mation  algorithms. 

We  continued  to  develop  our  ARMA  modeling  approach  to  the  multiple 
source  location  problem.  Some  of  the  key  theoretical  issues  related  to  unique- 
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ness  were  addressed  and  partially  resolved,  and  a  specific  algorithm  was 
developed  and  tested.  This  approach  is  a  first  attempt  to  handle  simultaneously 
several  sources,  rather  than  perform  single  source  detection  and  localization 
and  then  be  faced  with  the  association  problem. 

Several  advances  were  made  in  the  area  of  applying  image  reconstruction 
techniques  to  the  multiple  source  location  problem.  An  efficient  minimum  vari¬ 
ance  estimator  was  developed,  which  provides  a  key  step  towards  a  practical 
implementation  of  such  algorithms  in  a  DSN  environment.  Solutions  were 
obtained  for  the  problem  of  joint  estimation  of  source  distribution  and  attenua¬ 
tion.  This  problem  is  centred  to  the  application  of  reconstruction  techniques  to 
multi-source  tracking  by  acoustic  arrays. 


Intimation  Algorithms  far  DSN 

A  major  achievement  during  the  last  reporting  period  was  the  development 
of  ladder  forms  for  time-series  modeling  and  estimation.  The  simplicity  and 
efficiency  of  ladder  algorithms,  which  provide  an  exact  least  squares  solution  to 
the  multichannel  autoregressive  modeling  problem,  is  a  significant  break¬ 
through  in  the  areas  of  estimation  and  adaptive  signal  processing.  The  theory 
and  application  of  these  ladder  forms  is  expected  to  have  a  strong  impact  on  the 
general  field  of  signal  processing  and  on  the  acoustic  signal  processing  of  the 
DSN  type,  in  particular.  A  tremendous  amount  of  interest  was  generated  by  our 
presentations  of  the  recursive  ladder  forms,  and  the  number  of  publications  in 
this  area  seem  to  be  increasing  at  an  exponential  rate.  A  byproduct  of  this  was 
the  development  of  a  sensitive  event  detection  technique  for  DSN  applications, 
and  its  demonstration  in  a  speech  processing  problem  (in  the  context  of  pitch 
detection). 

Another  very  significant  achievement  was  the  development  of  doubling  algo- 
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rithms  for  the  efficient  inversion  of  Toeplitz  matrices.  The  need  for  inverting 
Toeplitz  matrices  arises  in  numerous  problem  areas.  Our  technique  provides  a 
significant  reduction  in  the  amount  of  computation  needed  to  invert  large 
matrices  of  this  type.  This  development  will  have  an  impact  on  adaptive  array 
processing  and  adaptive  beamforming,  for  example.  Both  of  these  problems  are 
of  great  importance  to  DSN  systems.  Other  adaptive  signal  processing  tech¬ 
niques  that  are  needed  in  the  uncertain  environment  in  which  the  DSN  operates 
were  developed  as  part  of  our  research  effort. 

Development  Of  ASgnal  Processing  Workbench 

An  essential  part  of  the  process  of  developing  new  techniques  for  source 
parameter  estimation  is  the  testing  of  proposed  algorithms  by  computer  simula¬ 
tion.  Simulations  allow  us  to  study  various  issues  related  to  the  performance  of 
the  algorithms  and  to  gain  insights  into  their  behavior.  Analysis  alone  is  not 
adequate,  especially  when  recursive  stochastic  algorithms  are  concerned. 
Simulations  are  an  invaluable  tool  for  studying  issues  such  as  convergence  rates, 
robustness,  and  estimation  accuracy.  Therefore,  establishment  of  a  signal  pro¬ 
cessing  workbench  that  will  allow  easy  comparison  and  modification  of  process¬ 
ing  steps  and  data  sequences  is  important.  In  particular,  this  will  facilitate  the 
exchange  of  processing  tools  and  data  bases  among  DSN  researchers  and  max¬ 
imize  the  impact  of  our  research  results  on  the  design  of  a  DSN  system. 

In  the  proposed  workbench,  the  processing  modules  will  function  in  a  stand 
alone  fashion  with  standard  input /output  formats  such  that  complex  signal  pro¬ 
cessing  functions  are  built  by  the  simple  cascading  of  modules.  In  the  UNDC 
operating  system,  precompiled  stand  alone  modules  can  be  connected  together 

using  a  system  command  called  a  pipe.  The  pipe  connects  the  output  of  one 
^  Sm  Bull  SysUm  ftctoncoi  Journal,  YoL  57,  No.  8,  July-August  1978. 
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tnodule  to  the  input  of  the  next.  Since  each  module  is  a  precompiled  program, 
the  signal  processing  can  be  changed  quickly  without  recompiling  or  linking,  just 
by  forking  to  the  appropriate  module.  Using  the  UNIX  system  debugger  com¬ 
mands  allows  the  internal  parameters  of  a  signal  processing  module  to  be 
observed  and  modified  as  data  is  being  processed.  This  facility  to  monitor  and 
change  signal  processing  parameters  while  data  is  being  processed  enhances  the 
debugging,  development,  and  understanding  of  new  algorithms.  The  ability  to 
control  multiple  processes  is  not  available  in  the  debuggers  provided  with  the 
current  releases  of  the  UNIX  operating  system.  Some  experience  with  the  prob- 
em  of  adding  this  capability  has  been  obtained  on  the  DEC-1 1/34  through  the 
development  of  our  own  debugger.  This  experience  will  be  applied  to  extend  the 
available  debuggers  on  the  VAX.  We  are  continuing  to  develop  a  signal  process¬ 
ing  workbench  based  on  the  above  ideas. 

In  addition  to  these  main  thrust  areas,  progress  was  made  in  our  basic 
-esearch  into  distributed  alg.  -ithms.  Interesting  results  were  obtained  on  two 
limensional  systems,  which  have  implications  for  planar  sensor  arrays  of  the 
ype  encountered  in  a  DSN  system.  New  efficient  estimation  algorithms  were 
[eveloped  for  nearly  stationary  processes.  This  development  opens  the  way  for 
elaxing  the  stationarity  assumption  that  is  commonly  made  in  current  techr 
iques  for  location  estimation  and  spectral  estimation.  The  framework  related 
>  nearly  stationary  processes  lends  Itself  naturally  to  distributed  implementa- 
on  and  is,  therefore,  well  suited  to  DSN  applications. 

In  the  next  phase  of  the  project  we  plan  to  continue  our  investigation  on  two 
vela:  basic  research  into  fundamental  mathematical  questions  related  to  DSN 
stem  design,  and  the  development,  coding,  and  testing  of  signal  processing 
odules  for  a  prototype  DSN  system.  An  important  part  of  our  work  will  be  the 
sting  of  new  signal  processing  techniques  on  data  provided  by  Lincoln  Lab  ora- 
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tories.  This  will  enable  us  to  demonstrate  the  performance  gains  that  can  be 
achieved  by  more  sophisticated  processing  of  multi-sensor  data  in  a  DSN 
scenario. 

The  results  of  our  research  effort  were  summarized  in  numerous  publica¬ 
tions  as  indicated  by  the  list  in  Section  7.  The  objective  of  this  report  is  to 
present  an  overview  of  our  work  in  a  coherent  framework,  while  emphasizing  its 
significance  to  different  aspects  of  the  DSN  problem.  We,  therefore,  limit  our 
technical  discussions  to  brief  descriptions  of  the  main  results,  leaving  most  of 
the  detaUs  to  other  publications.  It  is  hoped  that  this  report  will  provide  an 
introduction  and  summary  of  the  research  performed  at  the  1SL  for  the  DSN 
project,  enabling  the  reader  to  identify  specific  topics  of  interest  and  pointing  to 
the  relevant  references  for  more  details.  For  completeness  we  have  included 
some  of  the  key  publications  as  appendices. 
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2.  SOURCE  LOCATION  TECHNIQUES 

L  Overview 

A  central  problem  in  a  distributed  sensor  network  is  the  determination  of 
e  presence  and  location  of  sources  of  acoustic  and  electromagnetic  energy  in 
e  area  under  surveillance.  Multiple  sensors  provide  information  which  need  to 
i  processed  and  integrated  in  order  to  obtain  a  picture  of  the  outside  world, 
me  of  the  classical  approaches  to  this  problem  are:  high  resolution  spatial 
ectral  estimation  [CGK],  beamforming  for  sensor  arrays  [VT],  and  Kalman 
tering  techniques  and  various  multitarget  tracking  techniques  (see  e.g.,  the 
irvey  in  [YBS]).  Most  of  the  work  in  this  area  was  done  in  the  context  of  a  sin- 
e  sensor  (e.&..  a  radar  system)  or  a  group  of  sensors  acting  as  a  single  sensor 
.g.,  a  linear  antenna  array  or  a  hydrophone  array). 

The  DSN  problem  has  certain  features  which  do  not  seem  to  fit  well  these 
assical  approaches.  These  include:  the  geographical  dispersion  of  the  sensors, 
te  need  for  distributed  processing  (all  the  methods  mentioned  above  use  a 
jmpletely  centralized  processing  architecture),  the  presence  of  constraints  on 
te  amount  of  data  that  can  be  communicated  from  the  sensor  to  the  process* 
g  nodes  in  the  network,  and  the  large  number  of  sensors  and  sources  that  may 
i  involved.  In  view  of  these  differences  we  found  it  necessary  to  explore  new 
ichniques  for  solving  the  source  location  problem  that  seem  better  suited  for 
5N  applications,  including  bistatic  and  '’cooperative’'  modes. 

In  our  treatment  of  the  source  location  problem  we  considered  three  main 
sproaches.  The  first  Involves  the  estimation  of  time-differences-of-arrival 
DOA)  of  the  signal  to  several  sensors,  followed  by  estimating  the  location  of  the 
gnal  source.  This  type  of  approach  has  been  used  in  passive  sonar  systems, 
e  looked  much  more  closely  at  the  problem  structure  and  were  able  to  develop 
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some  improved  location  estimation  techniques  (Section  2.2  and  2.3).  Further¬ 
more,  this  analysis  led  to  a  characterization  of  source  identifiability,  i.e.,  under 

% 

what  conditions  is  it  possible  to  uniquely  locate  a  source  with  a  given  set  of 
measurements.  This  analysis  provides  some  suggestions  about  the  best  way  of 
placing  sensors  in  a  DSN  system. 

In  Section  2.4  we  show  that  the  source  location  problem  can  be  formulated 
as  a  system  identification  problem,  where  the  system  consists  of  several  sources 
and  sensors.  This  approach  provides  a  novel  way  of  simultaneously  estimating 
the  TDOA-s  and  spectral  characteristics  of  multiple  sources.  The  proposed  tech¬ 
nique  involves  multichannel  time-series  modeling.  Efficient  algorithms  for  per¬ 
forming  such  modeling  were  also  developed,  as  will  be  discussed  in  Section  3. 
These  algorithms  are  well-suited  for  distributed  implementation  and  provide  an 
attractive  alternative  to  FFT  related  techniques. 

In  Section  2.5  we  present  a  noncoherent  signal  processing  approach  for 
generating  an  image  of  the  area,  depicting  the  multiple  sources  that  are  present 
in  it.  This  approach  is  based  on  the  idea  of  image  reconstruction  from  line 
integral  projections,  which  has  been  used  extensively  in  medical  applications 
such  as  Computer  Aided  Tomography  (CAT).  We  show  that  by  treating  the  meas¬ 
urements  of  certain  types  of  sensors  as  line  integrals  through  an  "energy  inten¬ 
sity  map",  it  is  possible  to  reconstruct  an  image  showing  the  location  and  inten¬ 
sity  of  various  sources  of  acoustic /electromagnetic  energy  in  a  given  area. 

Finally,  in  Section  2.6  we  discuss  a  technique  for  remote  sensing  of  the 
atmosphere,  to  estimate  propagation  parameters.  Knowledge  of  such  parame¬ 
ters  can  be  used  to  improve  the  accuracy  of  the  location  estimates. 
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virement  geometry  commonly  employed  in  CT  applications  can  lead  to  a 
dramatic  reduction  in  the  total  computation  required  for  at  least  one  of  these 
methods,  the  minimum  variance  estimator,  thereby  making  it  more  competitive 
with  the  various  approximate  inversion  techniques  [WM].  In  DSN  applications, 
computational  efficiency  is  crucial,  because  of  the  large  amounts  of  data 
involved  and  the  temporal  resolution  required. 

Summary  of  [NUN] 

As  mentioned  earlier,  the  acoustic  imaging  problem  involves  estimation  of 
both  the  image  of  sources  and  the  attenuation  of  the  signals  from  the  sources  to 
the  sensors.  This  problem  was  treated  by  Nunes  [NUN]  in  the  context  of  medical 
imaging.  The  following  is  a  summary  of  the  topics  discussed  in  Nunes’  thesis. 

A  usual  approach  to  the  reconstruction  problem  in  computerized  tomogra¬ 
phy  consists  of  representing  the  image  as  a  finite  vector  and  the  measurements 
as  a  linear  operation  on  this  vector  plus  an  additive  noise.  This  modeling 
approach  has  several  advantages  over  Radon  transform-based  techniques  in  the 
sense  that  the  measurement  geometry  is  not  critical,  allowing  good  results  when 
evenly  spaced  measurement  profiles  cannot  be  obtained;  in  addition,  the  meas¬ 
urement  noise  is  part  of  the  model,  leading  to  the  use  of  estimation  techniques 
for  designing  reconstruction  algorithms. 

We  considered  three  topics  associated  with  this  approach.  First,  the  model¬ 
ing  of  the  measurements  is  studied  when  the  elements  of  the  image  vector  are 
the  first  n  coefficients  of  the  expansion  of  the  three  dimensional  distribution  to 
be  reconstructed  into  a  complete  set  of  orthogonal  functions.  Then  several  algo¬ 
rithms  that  are  currently  used  for  transmission  reconstruction  are  studied  and 
new  results  are  presented.  Finally,  the  more  complex  problem  of  reconstruct¬ 
ing  an  emission  distribution  embedded  in  a  medium  with  an  unknown  absorption 


it  first  we  present  a  brief  summary  of  the  relevant  research  results. 


immary  of  [WM] 

In  computerized  tomography  the  number  of  x-ray  photons  transmitted 
irough  an  object  of  interest  is  measured  by  an  array  of  detectors.  This  array  is 
stated  by  a  fixed  angular  increment  between  each  set  of  measurements.  From 
lese  measurements  the  total  attenuation  along  each  measurement  path  is 
erived  and  an  image  of  the  plane  of  the  measurement  paths  may  be  recon- 
xucted.  A  detailed  description  of  this  application  and  several  approaches  to 
stations  are  discussed  in  [MAC],  [GOR],  [NS]  and  [CBM]. 

Exact  deterministic  methods  for  the  reconstruction  of  an  image  of  the 
leasurement  plane  can  not  be  directly  applied  to  computerized  tomography  for 
ro  reasons.  First,  the  measurements  contain  the  inevitable  measurement 
oise  along  with  noise  due  to  counting  statistics.  In  addition,  the  measurements 
re  not  continuous,  but  are  taken  only  at  discrete  intervals  and  angular  posi- 
.ons.  Because  of  these  complications,  reconstructing  an  image  is  an  estima- 
.on.  and  not  an  inversion  problem.  Estimation  theory  results  have  been  used  in 
ais  application  by  Herman  [HL],  Rockmore  [ROC],  Wood  [WOO],  Fortes  [FOR], 
nd  Nunes  [NUN],  using  Bayesian  estimation,  maximum  likelihood  estimation, 
linimum  variance  estimation,  and  maximum  a  posteriori  estimators  respec- 
.vely.  These  methods  generally  produce  results  that  are  superior  to  the  ones 
btained  via  approximate  inversion  techniques  currently  employed,  especially  in 
ases  of  low  dosage  and  high  measurement  noise  [ROC],  [WOO].  There  is  growing 
iterest  in  quantitative  and  automatic  processing  of  CT  images,  and  we  contend 
lat  statistically  optimal  images  are  more  appropriate  for  this  application. 

The  estimators  discussed  above  are  not  generally  used  in  practice  because 
f  the  computational  complexity  and  storage  requirements.  However  the  meas- 
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Thus  J  many  radars  are  observing  a  given  volume  of  space,  their  measurements 
represent  surface  integrals  through  the  radar  image  of  that  volume.  A  convolve 
and  backproject  reconstruction  scheme  for  3-D  surface  projections  is  presented 
in  [RDF].  In  the  case  of  a  ranging  radar  that  uses  an  omnidirectional  antenna, 
the  projections  will  be  on  the  surface  of  a  sphere.  Techniques  for  reconstruction 
from  spherical  projections  can  be  similarly  derived. 

A  somewhat  different  formulation  has  to  be  used  in  the  case  of  acoustic  sen¬ 
sors.  An  acoustic  array  can  form  multiple  beams  in  different  spatial  directions. 
The  energy  received  at  the  array  output,  is  the  integrated  acoustic  energy  in  the 
beam  volume.  Thus,  a  system  consisting  of  several  acoustic  arrays  observing  a 
given  area,  provides  a  set  of  line  integral  measurements,  in  a  fan-beam 
geometry.  These  measurements  can  be  used  to  reconstruct  the  acoustic  energy 
map  of  the  observed  area,  from  which  the  source  locations  can  be  estimated. 
The  situation  is  somewhat  more  complicated  than  a  standard  reconstruction 
problem,  because  of  the  R2  attenuation  (/?  =  range)  due  to  spreading  loss. 

Several  aspects  of  the  image  reconstruction  approach  were  addressed  in 
our  research  program.  The  work  of  S.  Wood  [WM]  treated  the  computational 
aspects  of  image  reconstruction  algorithms,  and  derived  a  fast  implementation 
of  the  minimum  variance  estimator.  The  work  by  Nunes  [NUN]  and  Fortes  [FOR] 
addressed  the  joint  estimation  of  emission  and  attenuation  characteristics. 
These  type  of  solutions  can  be  used  to  locate  sources  in  the  presence  of  various 
types  of  propagation  loss.  Host  of  the  work  on  image  reconstruction  deals  with  a 
two-dimensional  situation.  This  work  also  extends  the  analysis  to  three  dimen¬ 
sions.  making  it  possible  to  address  the  problem  of  locating  sources  both  in 
range /bearing  and  in  altitude.  The  research  mentioned  above  was  done  in  the 
context  of  medical  applications.  However,  this  approach  is  very  promising  for 
solving  multiple  source  location  problems  a s  discussed  in  the  following  section. 


-  18- 


25  Image  Reconstruction  Techniques  For  Multiple  Source  Location 


2.5.1  Introduction  and  Summary 

A  class  of  inverse  problems  which  appears  in  many  scientific  and  engineer¬ 
ing  applications  is  the  reconstruction  of  a  density  function  from  projections.  In 
Computer  Assisted  Tomography  (CAT),  x-rays  or  7-rays  are  used  for  imaging  the 
internal  structure  of  the  body.  In  radio-astronomy,  the  electromagnetic  energy 
measured  by  antenna  arrays  is  used  to  reconstruct  the  radio-brightness  map  of 
the  celestial  sphere.  Various  techniques  have  been  developed  for  solving  the 
image  reconstruction  problem  and  this  area  is  by  now  fairly  well  developed  [BD] 
[SCU]  [GH].  In  this  section  we  will  briefly  indicate  how  image  reconstruction 
techniques  can  be  applied  to  the  DSN  problem.  See  also  [FDR]  for  a  more 
detailed  discussion. 


Let  us  consider,  for  example,  a  radar  system  observing  isotropic  point 
scatterers  distributed  in  a  volume  of  space.  Typically,  its  range  resolution  will 
be  much  better  than  its  azimuth  or  elevation  resolution.  The  radar  will  measure 


the  total  energy  of  returns  from  a  volume  of  space  which  is  essentially  a  thin 


Figure  22  Surface  Projection  in  Radar  Measurements 
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original  presentation  of  this  approach  we  suggested  the  multichannel  Recursive 
Maximum  likelihood  technique.  However,  the  development  of  computationally 
efficient  and  robust  ladder  algorithms,  leads  us  to  prefer  them  for  actual  imple¬ 
mentation.  Further  motivation  for  this  choice  is  some  recent  work  which  indi¬ 
cates  that  a  version  of  the  covariance  ladder  form  can  be  used  to  go  from  Left- 
to-Right-MFD,  a  crucial  step  in  our  tracking  technique.  Furthermore,  the  ladder 
algorithms  seem  to  be  well  suited  for  distributed  implementation,  and  thus  for 
application  in  a  DSN  system. 


A  .) 


‘j»  ^  /  ■> 


-  16- 

earthquake  localisation,  intruder  detection  or  artillery  localization).  The  active 
tracking  problem  (radar,  active  sonar)  can  also  be  handled  in  this  framework. 

The  proposed  approach  is  based  on  the  idea  of  fitting  a  multi-input  multi¬ 
output  model  to  the  vector  time-series  observed  at  the  outputs  of  the  sensors. 
Under  certain  assumptions  the  parameters  of  this  model  can  be  shown  to  con¬ 
tain  information  about  the  locations  of  all  the  sources,  as  well  as  other  informa¬ 
tion  such  as  source  spectra.  If  the  parameter  estimation  is  performed  recur¬ 
sively  and  continually,  it  is  hoped  that  the  relationship  between  the  model 
parameters  and  the  sources  to  which  they  correspond  will  be  consistently  main¬ 
tained.  Once  the  sources  are  labeled  in  a  particular  way,  this  labeling  will  stay 
fixed  without  need  for  rechecking  or  relabeling.  This  will  eliminate  the  need  for 
a  separate  step  of  association,  which  is  inherent  in  current  multi-source  track¬ 
ing  techniques. 

The  formulation  of  the  tracking  problem  as  a  multichannel  AKMA  modeling 
problem  was  presented  in  [FMl]  and  will  not  be  repeated  here.  We  note,  how¬ 
ever,  that  the  formulation  of  source  tracking  as  a  multichannel  signal  modeling 
problem  raises  a  number  of  interesting  questions  in  the  areas  of  system 
identification  and  structure  of  multivariable  linear  systems.  These  questions 
include: 

How  to  find  a  unique  (not  necessarily  minimum  phase!)  spectral  factor, 
under  structural  constraints; 

How  to  estimate  directly  the  parameters  of  a  right  Matrix  Fraction  Descrip¬ 
tion  from  a  given  data  set; 

How  to  include  certain  structural  constraints  into  the  multichannel  system 
identification  problem. 


Various  algorithms  can  be  used  to  estimate  the  ARMA  parameters.  In  our 
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2.4  An  AKICA  Modeling  Approach  to  Locating  Multiple  Sources  [Fill] 

fracking  multiple  sources  represents  special  difficulties  since  there  can  be 
uncertainties  associated  with,  the  measurements  beyond  their  inaccuracy,  which 
is  usually  modeled  by  some  additive  noise.  This  additional  uncertainty  is  related 
to  the  origin  of  the  measurements.  Since  several  sources  are  present,  it  is 
necessary  to  sort  out  which  measurement  corresponds  to  which  source.  In 
other  words,  in  addition  to  the  problem  of  detection  and  bearing /range  estima¬ 
tion.  there  is  a  problem  of  properly  labeling  the  set  of  measurements.  The 
latter  problem  is  usually  referred  to  as  source  association  or  track  formation. 

Typically,  these  two  facets  of  multi-source  tracking  are  treated  separately. 
First,  a  set  of  potential  source  locations  is  obtained.  Then  some  method  is  used 
to  label  these  locations  by  the  sources  to  which  they  correspond  in  a  manner 
consistent  with  previous  measurements.  Techniques  for  labeling  or  multi-source 
tracking  have  been  developed  using  various  approaches  including:  Kalman 
filtering,  Bayesian  methods,  Integer  programming,  and  Track  splitting. 

In  all  of  these  techniques  the  basic  detection  and  location  estimation  are 
performed  separately  for  each  source.  The  multi-source  aspect  of  the  problem 
enters  only  through  the  labeling  procedure.  In  other  words,  the  bracking  prob¬ 
lem  is  treated  as  a  collection  of  single  source  problems,  which  has  to  be  put 
together  in  a  systematic  and  consistent  way. 

In  [FMl],  we  attempt  to  tackle  directly  the  multichannel  nature  of  the  prob¬ 
lem.  Instead  of  looking  at  one  source  at  a  time,  we  want  to  estimate  simultane¬ 
ously  multi-source,  parameters  (such  as  location  and  spectrum).  The  approach 
is  beet  understood  in  the  context  of  the  passive  tracking,  where  an  array  of  sen¬ 
sors  measures  signals  (electromagnetic,  acoustic  or  seismic)  generated  by 
sources.  This  type  of  problem  arises  in  sonar  systems,  acoustic  surveillance  sys¬ 
tems  (detection  of  low  flying  aircraft)  and  seismic  arrays  (oil  exploration, 
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Gaussian  second  order  filter  developed  by  Sorenson  and  Stubbard  [SS]  using 
some  additional  assumptions  on  the  structure  of  the  measurement  update  equa¬ 
tions.  Our  derivation  is  more  general  in  that  it  does  not  make  any  a  priori 
assumptions  regarding  the  structure  of  the  estimator.  The  second  order  Gaus¬ 
sian  filter  is  known  to  have  better  convergence  properties  than  the  EKF  and  is, 
therefore,  preferable  for  tracking  applications. 

The  second  order  Gaussian  filter  lends  itself  to  distributed  implementation. 
Several  distributed  forms  of  the  filter  equations  are  derived  in  [LM3].  A  com¬ 
parison  with  a  completely  centralized  scheme  demonstrates  that  our  distributed 
algorithms  lead  to  considerable  savings  in  communication  requirements.  More¬ 
over,  they  exhibit  a  parallel  structure  and  are  flexible.  We  believe  that  these 
algorithms  will  provide  a  practical  intermediate  solution  for  tracking  in  a  DSN 
system.  However,  a  solution  based  on  the  conditionally  linear  formulation  ear¬ 
lier  should  be  considered  as  the  ultimate  goal. 
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**♦1  s  Fk*k  *  2  e  **,;$**  +  **  (2.3) 

i*i 

where  $el,  e*  ....  *”}  is  the  standard  basis  of  Rn  and  the  matrices  C&  are 
symmetric. 

The  presence  of  a  quadratic  measurement  equation  leads  to  a  nonlinear 
estimation  technique  which  is  typically  solved  by  using  an  Extended  Kalman 
Filter  (EKF)  [JAZ].  The  EKF  is  based  on  a  successive  linearization  of  the  meas¬ 
urement  and  state  equations.  While  the  EKF  seems  to  work  well  in  some  track¬ 
ing  situations,  various  authors  report  encountering  serious  difficulties  especially 
in  near  range  and  low  signal-to- noise  situations.  A  principal  drawback  of  the  EKF 
is  that  it  is  not  guaranteed  to  converge.  Divergence  of  the  estimation  error  is 
fairly  common,  especially  under  the  conditions  mentioned  above. 

The  problem  formulation  we  presented  in  Section  2.2  and  in  [LM4],  opens 
the  way  for  alleviating  these  difficulties.  In  our  approach  the  measurement 
equations  are  linear  (or  at  least  conditionally  linear).  If  in  addition,  the  source 
dynamics  are  assumed  to  be  also  linear  (a  common  assumption),  the  problem 
can  be  solved  using  the  standard  Kalman  Alter!  This  circumvents  the  need  for 
any  linearizations  and  the  associated  inaccuracies  and  divergence  problems. 
The  combination  of  our  static  source  location  technique  with  a  linear  Kalman 
Alter  has  the  potential  of  providing  an  elegant  solution  to  the  tracking  problem. 
We  should  note,  however,  that  this  method  is  not  yet  fully  developed  and  we  are 
still  working  on  some  of  the  issues  that  need  to  be  resolved  before  a  practiced 
application  can  be  attempted. 

Another  way  of  alleviating  divergence  problems  is  to  develop  an  alternative 
nonlinear  estimation  technique.  Under  the  assumption  of  a  Gaussian  a-posterior 
density  function,  we  obtained  in  [LM3]  a  finite  order  solution  that  does  not 
involve  linearization.  Our  estimation  algorithm  turns  out  to  be  identical  with  the 
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2.3  Distributed  Estimators  Ftor  Quadratic  Systems  [LM3] 

The  estimation  technique  described  in  the  previous  section  was  presented 
as  the  solution  to  a  static  problem:  given  a  set  of  TDOA  measurements  * 

corresponding  to  a  certain  source  location,  find  the  best  estimate  of  that  loca¬ 
tion.  The  source  tracking  problem  involves,  however,  a  dynamic  aspect  that  is 
not  captured  by  the  discussion  presented  earlier.  In  this  section  we  present  an 
approach  that  addresses  the  time-varying  features  of  tracking  one  or  more 
sources. 

A  standard  approach  to  the  tracking  problem  is  to  construct  a  nonlinear 
dynamic  equation  that  describes  the  evolution  of  the  source  state  vector.  This 
state  vector  will  typically  contain  the  source  position  coordinates  (two  or  three 
variables),  its  velocity  and  perhaps  the  clock  bias.  See  for  example,  the  discus¬ 
sions  in  [NM]  -  [DMFl].  The  measurement  equation  can  be  shown  to  exhibit  on 
intrinsic  quadratic  structure 

+  vh  ,  (2.1) 

where 

I  =  (1.  1 . lj  x  RJ> 

=  white  noise  process 

The  measurement  matrix  Hk  may  be  time-varying,  when  the  sensors  them¬ 
selves  are  moving.  The  state  vector  typically  obeys  a  linear  evolution  equation 

of  the  type 

**♦1  =  Fkxh  +  ti*  .  (2.2) 

where  «*  is  a  white  noise  process,  and  Fk  is  a  matrix  representing  source 
motion  dynamics.  In  fact,  the  approach  can  handle  with  no  significant  increase 
of  complexity  nonlinear  dynamics  of  a  quadratic  type,  l.e., 
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[LM3]. 

The  use  of  a  two  step  procedure,  i.e.,  involving  a  preprocessing  step,  is  very 
general  for  this  estimation  problem.  In  [SCH]  a  TDOA  "averaging"  procedure  was 
proposed  and  a  straightforward  generalization  of  the  idea  is  to  estimate  from 
the  “raw"  or  measured  TDOA’s  via  linear  least  squares  techniques  the  ranges  up 
to  an  arbitrary  constant,  that  can  be  viewed  as  a  clock  bias;  a  consistent  set  of 
TDOA's  is  directly  obtained  from  these  estimates.  Since  more  structural  infor¬ 
mation  is  used  in  our  preprocessing  step  than  in  a  simple  TDOA  averaging  step, 
the  corresponding  estimates  eure  in  general  more  accurate.  The  overall  estima¬ 
tion  procedure  also  performs  in  general  better  than  more  traditional 
approaches,  especially  since  the  equation  coefficients  are  measurement 
independent. 

In  [DMFl]  the  convenient  approach,  which  separates  the  TDOA  estimation 
and  the  source  location  problems  has  been  followed.  However,  it  is  clear  that  a 
completely  satisfactory  answer  to  both  the  identiflability  and  estimation  prob¬ 
lems  will  only  result  from  a  global  approach.  A  close  combination  of  the 
approaches  of  [DUFl]  [DMF2]  and  [FM1]  appears  to  be  possible  and  we  are  con¬ 
tinuing  our  investigation  in  that  direction 


[mi]  *nd  [m]. 


In  [DMF2]  we  first  determine,  as  described  above,  the  joint  source  and  sta¬ 
tion  geometries  for  which  ambiguities  in  the  determination  of  the  source  loca¬ 
tion  arise.  Then,  using  this  result,  we  are  able  to  introduce  a  simple  rule  for  the 
placement  of  the  stations  to  ensure  that  the  source  coordinates  are  always 
identifiable.  It  turns  out  that  a  linear  array  does  not  satisfy  that  rule  and  indeed 
the  sensitivity  of  any  source  location  estimate  (geometric  dilution  of  precision) 
is  very  high  when  the  source  lies  close  to  the  axis  of  the  array.  In  the  remainder 
of  the  paper  the  geometry  of  the  stations  is  assumed  to  satisfy  the  placement 
rule.  Under  this  mild  condition  our.  pxp  matrix  has  additional  properties  that 
allow  us  to  devise  several  simple  estimation  procedures. 

Our  estimation  procedures  involve  two  distinct  steps  and,  assuming  infor¬ 
mation  from  the  second  step  is  not  fed  back  into  the  first  step,  are  therefore 
subop timal;  their  advantage  lie  in  their  relative  simplicity.  The  first  step,  best 
viewed  as  a  preprocessing  step,  consists  of  the  (linear  least  squares)  estimation 
of  the  ranges  between  the  source  and  the  stations.  In  [DMFl]  a  simple  function 
of  the  ranges  was  estimated  rather  than  the  ranges  themselves.  This  estimate 
was  more  robust  but  led  to  a  more  complex  second  step.  Note  that,  at  this  level, 
the  range  estimates  could  be  improved  if  some  direct  range  measurements  sure 
available.  In  the  second  step  the  pxp  matrix  is  factored  using  this  range  infor¬ 
mation;  this  yields  a  set  of  linear  and  quadratic  equations  in  the  source  coordi¬ 
nates  that  can  be  viewed  as  our  new  or  refined  observations,  replacing  the 
TDOA's.  The  coefficients  appearing  in  these  equations  are  not  a  function  of  noisy 
measurements,  they  are  functions  of  the  stations  coordinates  ordy.  Because 
this  property  is  satisfied,  standard  estimation  procedures  can  be  applied  to 
these  equations,  in  the  static  case  we  can  use  for  Instance  the  results  in  [HB],  or 
if  a  dynamic  model  is  available,  a  second-order  filter  can  be  used,  see  e.g., 


A  reasonable  way  to  resolve  such  a  difficulty  is  to  reduce  the  original  equa¬ 
tions  to  an  "equivalent"  set  of  equations  for  which  known  estimation  procedures 
are  more  easily  applied.  "Equivalent”  means  that  the  mapping  from  one  set  of 
equations  to  the  other  is  one-to-one,  i.e.,  no  information  is  lost  through  that 
mapping.  R.  Schmidt  made  an  interesting  attempt  in  that  direction  in  [SCH]. 
However,  we  have  shown  in  [DMF1]  that  the  linear  equations  he  obtained  are 
highly  redundant  and  contain  less  information  than  the  original  TDOA’s. 

Simple  notions  of  spatial  rotation  and  translation  invariance  led  us  to  con* 
aider  a  particular  pxp  matrix.  The  (ij)-th  entry  is  a  function  of  the  distance 
and  range  difference  between  the  (i.j)-th  pair  of  stations  and  has  a  simple 
expression  in  terms  of  the  stations  and  source  positions.  That  matrix  equation 
plus  some  auxiliary  sign  information  are  easily  shown  to  be  equivalent  to  the  ori¬ 
ginal  set  of  range  differences. 

An  immediate  advantage  of  rewriting  of  the  TDOA’s  equations  is  that  they 
can  be  analysed  systematically  using  the  simple  tools  of  linear  algebra.  The 
cases  for  which,  given  the  stations  position  and  (noise  free)  range  differences, 
the  equations  have  more  than  one  solution  are  completely  characterized  in  the 
2-  and  3-dimensional  Euclidean  spaces.  Thus  the  static  "deterministic” 
identiflability  problem  is  solved.  This  can  be  viewed  as  a  first  step  in  the  solution 
of  the  "stochastic"  identiflability  problem,  i.e.,  given  noisy  measurements,  for 
both  the  static  and  moving  source  cases.  Such  a  study  is  of  fundamental 
interest  since  it  indicates  how  much  information  about  the  source  location  can 
be  extracted  from  the  TDOA’s  and  therefore  the  limitations  of  an  estimation  pro¬ 
cedure.  We  carried  out  such  a  study  in  [DMF2].  The  stochastic  identiflability 
problem  1b  discussed  there  along  the  lines  of  [PIC].  Ways  to  extend  the 
approach  to  multi-source  and  multi-path  situations  are  also  Indicated;  however, 
the  detailed  study  of  these  cases  is  not  carried  out  in  this  paper  and  we  refer  to 


2.2  Source  Location  from  Time  Differences  of  Arrival  [DMF1]  [DMF2] 


The  problem  of  source  location  from  time  differences  of  arrival  (TDOA)  is  a 
standard  nonlinear  estimation  problem.  The  usual  approach  to  source  location 
estimation  is  based  on  the  construction  of  hyperbolic  lines  of  position  (LOP). 
Each  measurement  of  TDOA  at  a  pair  of  stations  (or  after  multiplication  by  the 
velocity  of  propagation,  each  range  difference)  determines  an  LOP  that  is  a 
branch  of  a  hyperbola.  Then  the  source  location  may  be  inferred  from  the  inter¬ 
sections  of  the  j^]  LOP’s  obtained  from  p  stations.  See  Figure  2.1  for  the  case 
of  3  stations. 


Figure  2.1.  Hyperbolic  lines  of  Poeitian 

Unfortunately,  information  about  the  measurement  noise  cannot  be  completely 
accounted  for  in  this  approach,  so  that  it  seems  impossible  to  claim  some 
optimality  for  any  estimator,  or  to  characterize  the  statistical  properties  of  the 
associated  estimate. 
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coefficient  distribution  is  treated  as  a  joint  estimation  problem  where  the  meas¬ 
urement  equation  is  no  longer  linear. 

The  first  part  of  this  thesis  develops  an  expression  for  transmission  meas¬ 
urements  that  takes  into  account  the  quantization  noise  and  the  mathematical 
approximations  necessary  for  the  linearization  of  the  measurement  equation. 
Equations  for  emission  and  absorption  measurements  are  also  derived  in  this 
part. 

Considering  transmission  reconstruction  techniques,  it  is  shown  that  two 
recently  proposed  algorithms  asymptotically  converge  to  the  weighted  least 
square  solution.  It  is  also  shown  that  a  priori  information  about  the  image  vec¬ 
tor  can  be  included  in  the  measurement  set  in  order  to  make  the  algorithms 
converge  to  the  linear  minimum  variance  solution.  An  additional  result  is  that 
the  well  known  Algebraic  Reconstruction  Technique  (ART)  converges  to  a 
weighted  least  squares  solution  if  a  suitable  time  varying  relaxation  coefficient  is 
used.  We  also  show  that  convergence  is  achieved  even  if  the  set  of  measure¬ 
ments  is  not  consistent. 

Regarding  emission  reconstruction  problems  in  the  presence  of  unknown 
absorption,  a  nonlinear  algorithm  for  jointly  estimating  the  two  distributions  is 
developed,  leading  to  the  establishment  of  an  estimation  criterion  which  takes 
into  account  a  priori  information  about  the  image.  Additional  algorithms  are 
proposed  and  it  is  shown  that  they  converge  to  a  solution  minimizing  this  new 
criterion. 

Computer  simulations  of  the  transmission  and  joint  estimation  algorithms 
are  presented,  permitting  the  evaluation  of  their  performance  as  a  function  of  a 
limited  number  of  iterations. 


Summary  of  [FOR] 

In  a  related  work,  Fortes  [FOR]  develops  an  estimation  approach  to  3-D 
reconstruction  problems  involving  counting  statistics,  with  applications  to  medi¬ 
cal  imaging.  The  following  is  a  summary  of  the  topics  discussed  in  Fortes' 
thesis. 

The  problem  of  estimating  three  dimensional  distributions  from  a  given  set 
of  measurements  that  are  noisy  functionals  of  these  distributions  is  analysed. 
The  solution  to  this  problem  depends  on  several  factors,  but  mainly  on  the 
measurement  model  and  geometry. 

In  many  applications  measurement  noise  has  a  counting  nature  that  allows 
a  discrete  stochastic  process  characterization.  Hence,  a  general  model  assum¬ 
ing  a  Doubly  Stochastic  Poisson  Process  characterization  for  the  measurement 
functions  is  presented.  Based  on  this  model,  an  estimation  approach  leading  to 
the  Maximum  A  Posteriori  probability  (MAP)  estimate  is  carried  out. 

The  presented  formulation  is  clearly  applicable  to  a  large  family  of  three 
dimensional  reconstruction  problem&  However,  to  illustrate  our  results,  the 
medical  imaging  problem  of  Computerized  Tomography  is  mainly  used.  Com¬ 
puter  simulations  show  that  in  transmission  tomography  our  approach  leads  to  a 
non-linear  estimator  that  represents  an  improvement  over  the  linear  minimum 
mean  squared  error  estimate.  In  particular,  the  squared  error  obtained  with 
the  MAP  estimate  is  shown  to  have  a  lower  mean  and  a  narrower  distribution 
than  the  (me  obtained  with  the  Kalman  Filter.  In  the  case  of  emission  tomogra¬ 
phy  with  unknown  attenuation,  simultaneous  estimation  of  source  and  attenua¬ 
tion  distributions  is  obtained  using  the  MAP  technique.  In  both  cases,  the  per¬ 
formance  of  the  MAP  estimator  is  studied  through  a  comparison  with  the 
Cramer-Rao  mean  squared  error  lower  bound. 
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The  necessary  conditions  for  uniqueness  of  the  MAP  estimate  for  transmis¬ 
sion  and  emission  tomography  with  known  attenuation  are  determined  and 
shown  to  be  related  to  the  non-singularity  of  the  projection  matrices. 

2.5.2  Image  Reconstruction  Techniques  Applied  to  Target  Detection 

In  Medical  Imaging  and  Radioastronomy  the  images  are  generally  "com¬ 
pact"  and  the  ratio  between  the  required  resolution  and  the  size  of  the  image  is 
not  excessive,  in  the  target  detection  problem  the  images  are  very  sparse  and 
the  resolution  required  very  large.  This  means  that  an  excessive  number  of 
basis  functions  will  be  necessary  in  order  to  discretize  the  image.  If  the  usual 
medical  image  reconstruction  techniques  are  applied  directly  to  the  target 
detection  problem,  a  colossal  problem  in  terms  of  memory,  computational  time 
and  number  of  required  measurements  will  result.  If  the  procedures  are 
modified  to  consider  the  sparseness  and  positiveness  of  the  images  then  medi¬ 
cal  imaging  reconstruction  techniques  become  efficient  and  valuable  tools  for 
solving  of  the  multi-target  detection  problem.  There  are  at  least  two  ways  of 
making  medical  Image  reconstruction  techniques  suitable  for  multi-target 
detection:  pre-processing  the  data  in  order  to  reduce  the  actual  image  to  a  set 
of  regions  which  are  likely  to  contain  targets,  or  dividing  the  image  in  pixels 
much  greater  than  the  expected  size  of  a  target  and  estimating  the  total  target 
intensity  for  each  pixel. 

In  medical  imaging  the  image  is  generally  divided  in  n  adjacent  nonover¬ 
lapping  pixels,  which  can  be  seen  as  the  first  n  basis  functions.  If  the  number  of 
pixels  is  large  enough,  the  quantization  noise  can  be  neglected.  Of  course  the 
target  detection  problem  can  be  modeled  the  same  way',  however,  considering 
that  the  size  of  a  target  is  very  small  compared  to  the  region  under  surveillance, 
a  very  large  number  of  pixels  will  be  required  in  order  to  keep  the  quantization 
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noise  down.  This  implies  a  large  amount  of  memory  and  computation  time.  For 
linear  reconstruction  techniques,  at  least  as  many  measurements  as  the 
number  of  pixels  should  be  used  in  order  to  have  a  unique  solution,  using  the 
least-squares  criterion  for  instance.  One  can  argue  that  the  linear  minimum 
variance  solution  will  always  be  unique.  On  the  other  hand,  one  can  not  expect 
good  results  when  the  number  of  measurements  is  small  because  essentially 
only  second  order  prior  statistics  are  used  to  complement  the  measurement 
set,  but  second  order  priors  cannot  represent  two  major  a  priori  characteristics 
of  the  image:  positiveness  and  sparseness. 

Of  course  optimal  nonlinear  estimation  criteria,  like  M.M.S.E.  or  M.A.P.,  can 
account  for  positiveness  and/or  sparseness.  Unfortunately,  if  the  pixel  represen¬ 
tation  is  kept,  the  amount  of  memory  and  computation  time  is  still  prohibitive 
at  the  moment  without  special  purpose  VLSI  hardware.  On  the  other  band,  the 
sparseness  and  positiveness  of  the  image  intuitively  suggest  the  use  of  nonlinear 
procedures  as  a  preliminary  processing  of  the  image.  This  pre-processing  would 
determine  which  pixels  can  be  assumed  empty  and  exclude  them  of  further  pro¬ 
cessing.  Linear  techniques  can  then  be  used  for  processing  the  remaining  pixels. 
The  pre-processing  can  also  provide  a  priori  statistics  for  the  linear  processing. 

A  data  pre-processing  procedure  has  been  simulated  that  reduces  the 
dimensionality  of  the  multi-target  detection  problem.  The  method  uses  the 
sparseness  and  positiveness  characteristics  of  the  image  for  establishing 
"apriori"  mean  and  covariance  tor  the  Anal  processor.  Since  the  image  is  posi¬ 
tive  and  sparse,  it  is  reasonable  to  assume  that  a  measurement  close  to  zero  is 
nonzero  only  because  there  is  noise  involved,  and  all  the  pixels  combined  to 
form  that  measurement  can  be  set  to  zero.  Whether  a  measurement  should  be 
considered  zero  or  not  is  decided  by  comparing  it  with  a  threshold  that  is  a 
function  of  the  noise  level  and  of  the  other  measurements.  If  a  given  pixel  con- 
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tributes  to  nonzero  measurements  in  all  profiles,  then  we  assign  to  it  a  mean 
and  a  variance  which  give  the  Kalman  filter  (processor)  a  great  flexibility  in 
determining  its  value.  If  a  pixel  contributes  to  nonzero  measurements  in  all  but 
one  profile,  then  we  set  to  it  a  smaller  mean  than  before  and  a  smaller  variance, 
Le.  we  "tell"  the  processor  that  those  pixels  should  be  small. 

The  detection  of  targets  using  pre-processing  as  described  above  has  been 
computer  simulated.  The  model  for  the  true  signal  point  sources  of  differing 
intensity  (referred  to  as  phantom)  and  the  geometry  is  illustrated  in  Figure  2.3. 
The  additive  noise  was  assumed  constant.  Figure  2.4  illustrates  the  lack  of  noise 
suppression  in  the  backprojected  reconstruction.  The  preprocessing  with 
different  thresholds  is  illustrated  in  Figure  2.5.  The  number  and  amplitude  of  the 
spurious  sources  increases  as  the  threshold  increases. 

Another  approach  to  the  target  detection  problem  is  to  redefine  our  objec¬ 
tives.  Suppose  that,  instead  of  trying  to  determine  the  exact  position  of  a  target, 
we  divide  the  region  under  observation  in  several  pixels,  each  of  them  much 
bigger  than  the  expected  size  of  a  target.  We  now  want  to  determine  the 
integral  of  the  image  in  each  pixel,  instead  of  the  image  itself.  Of  course  this  set 
of  parameters  give  us  information  about  the  distribution  of  targets  in  the  area. 
This  information  could  be  sufficient  in  many  cases.  In  this  case  a  measurement 
is  not  an  operation  over  the  parameter  space  Itself,  but  over  the  target  distribu¬ 
tion  image.  Of  course  it  is  possible  that  a  linear  operation  over  the  image  can  be 
also  represented  as  a  linear  operation  over  the  parameter  space.  Loosely 
speaking,  if  all  points  of  a  pixel  are  treated  in  the  same  manner,  they  can  be 
grouped  together  and  replaced  by  the  value  of  the  integral  over  the  pixel.  In 
practice,  it  is  not  possible  to  treat  all  the  points  in  a  pixel  in  the  same  way.  How¬ 
ever  in  many  cases  we  can  say  that  they  are  treated  approximately  in  the  same 
way,  which  yields  an  approximated  expression  for  a  measurement  in  terms  of 
the  parameter  space  only. 
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2.6  Optical  Probing  of  the  Atmosphere 

The  DSN  concept  is  not  limited  to  any  particular  class  of  sensors.  However, 
the  demonstration  system  designed  by  the  Lincoln  Lab  group  will  make  use  of 
small  acoustic  arrays  to  detect  and  track  low  flying  aircraft.  The  performance 
of  such  a  system  will  depend  among  other  things  on  the  properties  of  the  propa¬ 
gation  medium.  i.e.,  the  atmosphere.  Knowledge  of  various  parameters  such  as 
temperature  profile  and  wind  velocity  can  improve  location  accuracy,  especially 
at  longer  ranges.  Working  on  a  different  problem,  S.  Delateur  [DEL]  investigated 
the  possibility  of  estimating  various  atmospheric  variables  based  on  an  optical 
probing  system.  The  remote  sensing  techniques  developed  in  his  thesis  may  be 
applicable  to  an  advanced  DSN  system.  We  next  summarize  the  work  described 
in  [DEL]. 

When  an  optical  wave  travels  through  clear  air  turbulence,  inhomogeneties 
in  the  index  of  refraction  cause  intensity  variations  at  a  receiving  or  imaging 
plane.  These  effects  allow  remote  sensing  of  the  earth's  atmosphere.  Two  meas¬ 
urements  of  interest  are  path  profiles  of  turbulence  strength  and  wind  velocity. 

The  lack  of  path  selectivity  and  difficulties  in  signal  processing  (due  to  the 
non-stationarity  of  the  atmosphere)  have  characterized  methods  which  employ 
covariance  calculations  of  the  received  field  to  estimate  cross-path  wind  veloci¬ 
ties.  With  a  single-scatter,  phase-screen  model  for  clear  air,  an  analysis  of  line- 
of-sight  propagation  was  presented  using  Fourier  optics,  including  a  discussion 
of  the  assumed  nature  of  the  turbulence  spatial  spectrum  Specifications  are 
detailed  for  the  transmitter  and  receiver  apertures  and  the  processing  of  the 
received  signal  which  provide  estimates  of  the  turbulence  strength  and 
transverse  wind  speed  at  selected  points  along  the  path. 


An  Incoherent  optical  system  incorporating  these  principles  has  been  con¬ 
structed.  Results  from  experiments  to  monitor  the  movement  of  the 
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atmosphere  simultaneously  at  four  turbulent  scale  sizes  are  reviewed  with 
direct  comparison  between  optical  measurements  of  wind  speed  and  local 
anemometer  readings  included. 

Using  computer  simultations  of  the  received  signal  characteristics,  the  sys¬ 
tem  problems  concerning  velocity  measurement  errors  and  path  profile  selec¬ 
tivity  were  explored.  Final  comments  contain  improved  signal  processing  tech¬ 
niques  and  aperture  specifications  for  future  equipment. 


3.  ESTIMATION  ALGORITHMS  FOR  DSN 


3.1  Overview 

The  problem  of  locating  and  classifying  multiple  sources  on  the  basis  of 
information  received  by  different  sensors  in  the  network  requires  the  use  of 
various  signal  processing  modules.  Many  of  these  modules  embody  some  kind  of 
an  estimation  algorithm.  The  performance  of  the  DSN  system  will  depend 
strongly  on  the  quality  of  its  "front  end",  i.e..  the  part  of  the  system  that  takes 
raw  sensor  information  and  converts  it  into  higher  level  parameters  such  as 
TDQA  estimates,  spectral  estimates,  or  likelihood  functions.  Therefore,  the 
development  of  estimation  algorithms  that  are  efficient  and  suitable  for  distri¬ 
buted  implementation  is  an  essential  step  towards  the  achievement  of  a  high 
performance  DSN  system.  Since  multiple  sensors  and  sensor  arrays  are 
involved,  we  need  in  general  to  consider  multichannel  estimation  algorithms. 

A  major  step  in  this  direction  is  provided  by  our  development  of  ladder 
forms  for  time-series  modeling  and  estimation.  This  development  can  be  con¬ 
sidered  as  a  significant  breakthrough  in  the  areas  of  least-squares  estimation 
and  adaptive  signal  processing.  The  body  of  theory  and  the  numerous  algo¬ 
rithms  developed  during  the  last  reporting  period  will  have  a  strong  impact  in 
many  areas,  including:  speech  processing,  acoustic  signal  processing  for  sonar, 
adaptive  arrays,  and  high  resolution  spectral  estimation.  A  tremendous  amount 
of  interest  was  generated  by  our  presentations  of  the  theory  and  practice  of 
Ladder  forms.  The  number  of  publications  and  applications  related  to  ladder 
forms  seems  to  be  on  an  exponential  growth  curve.  The  highlights  of  these  very 
exciting  results  are  presented  in  Section  3.2.  VTe  believe  that  estimation  tech¬ 
niques  based  on  these  ladder  forms  will  lead  to  significant  improvements  in  DSN 
performance,  over  that  achievable  with  the  more  traditional  techniques. 
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An  interesting  byproduct  of  the  ladder  form  was  the  fact  that  one  of  the 
quantities  being  ‘computed  as  part  of  their  recursions  is  a  sensitive  indicator  of 
events  in  the  observed  data.  This  led  us  to  develop  efficient  algorithms  that  can 
be  used  to  detect  the  appearance  or  disappearance  of  sources,  or  other  events 
in  a  DSN  scenario.  In  Section  3.3  we  summarize  some  of  the  results  in  this  area. 

Another  very  important  development  in  the  area  of  efficient  computational 
techniques  was  our  work  on  doubling  and  tree  algorithms.  One  main  result  of 
this  work  was  the  development  of  an  inversion  algorithm  for  NxN  Toeplitz 
matrices  that  requires  0(NlogzN)  operations.  The  need  for  inverting  Toeplitz 
matrices  arises  in  most  estimation  problems  related  to  stationary,  or  nearly  sta¬ 
tionary  (finite  rank)  processes.  For  general  matrices  we  developed  an  efficient 
distributed  algorithm,  that  requires  only  local  communications  between  proces¬ 
sors  either  arranged  in  a  one  or  two  dimensional  array.  An  application  of  partic¬ 
ular  interest  to  the  DSN  problem  is  adaptive  beamforming.  Section  3.4  presents 
a  brief  description  of  our  work  in  this  area. 

Ladder  forms  and  other  least-squares  estimation  techniques  have  numerous 
applications  to  adaptive  signal  processing  problems,  lie  investigated  several 
such  pro  We  ms  including  adaptive  line  enhancement  and  adaptive  linear  phase 
filtering.  Section  3.5  presents  these  and  other  related  results. 
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3.2  Ladder  Algorithms  For  Estimation  and  Adaptive  Processing 

A  large  number  of  publications  were  written  on  the  theory  and  application 
of  ladder  forms  during  the  last  reporting  period,  e.g.,  [LEE]  [LFM]  [LM1]  [FR2] 
[FR3]  [PFM].  The  number  and  diversity  of  research  topics  that  were  addressed 
makes  it  difficult  to  present  a  comprehensive  summary.  In  this  section  we  will 
therefore  discuss  only  the  highlights  of  our  work.  We  start  by  introducing  ladder 
structures  and  their  properties. 

Ladder  Structures 

Finite  dimensional  linear  systems  can  be  realized  using  different  network 
configurations.  These  are  typically  classified  as  direct,  cascade  and  parallel 
forms  [OS].  As  an  example,  consider  an  auto-regressive  moving-average  (ARMA) 
model  given  by 

VT  =  -  £  +  £  6*ur_,  .  (3.1) 

i«l  i*l 

A  direction  realization  of  this  model  is  depicted  in  Figure  3. 1. 


figure  3.1.  Direct  Realisation  of  H(z)  =  B[z)/A(z) 

Many  other  direct  realizations  of  the  same  transfer  functions  are  possible, 
corresponding  to  different  canonical  forms  (controller,  observer,  etc.  [KA1]). 
Cascade  realizations  have  been  used  extensively  in  network  synthesis  and 
related  areas  [W],  Perusal  of  the  vast  literature  on  adaptive  systems  reveals 
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t  an  overwhelming  majority  of  the  work  in  this  area  is  based  on  the  direct 
lization.  A  notable  exception  is  a  special  case  of  cascade  realizations:  the 
ier  form.  A  typical  example  of  a  ladder  form  for  an  autoregressive  (AR)  or 
pole  model  is  depicted  in  Figure  3.2. 


Figure  3.2.  An  Autoregressive  Ladder  Form 


der  type  structures  originated  in  Network  and  Scattering  Theory  and 
ame  popular  more  recently  in  speech  processing  [MG]  and  actually  earlier  in 
mic  applications.  However,  their  overall  impact  on  the  area  of  adaptive  pro¬ 
sing  is  still  relatively  small  when  compared  to  the  popularity  of  direct  reali- 

003. 

Ladder  structures  have  many  attractive  features  which  make  them  prefer- 
i  to  direct  realizations.  One  such  feature  is  the  orthogonality  (decoupling) 
perty:  the  signals  propagating  in  a  ladder  form  whitening  filter  (the  forward 

backward  innovations)  are  uncorrelated  and  therefore  the  coefficients  of 
filter  can  be  adjusted  independently.  One  manifestation  of  this  property  is 
fact  that  when  the  order  of  the  ladder  filter  is  increased,  one  can  add  addi- 
al  sections  to  the  filter  without  changing  the  rest  of  the  filter.  In  other 
is,  the  (p+l)-th  order  filter  which  fits  a  given  data  set  h/(0!<=o  .  is  the 
e  as  the  p-th  order  filter,  except  for  the  last  section,  i.e.,  the  lower  order 
iels  are  nested  in  the  higher  order  models.  Therefore,  ladder  filters  make  it 
lible  to  get  simultaneous  estimates  of  models  of  different  orders  (up  to  a 
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maximum  filter  order),  with  no  extra  computation!  This  flexibility  is  extremely 
useful  in  tackling  the  problem  of  model  order  determination,  which  is  one  of  the 
most  difficult  aspects  of  signal  modeling.  It  should  be  noted  that  this  nesting 
property  does  not  hold  for  direct  realizations. 

The  cascaded  structure  of  the  ladder  filter  (consisting  of  identical  sections) 
is  very  convenient  from  an  implementation  standpoint,  especially  if  special  pur¬ 
pose  hardware  (e.g.,  VLSI  )  is  used.  Other  properties  of  ladder  forms  which  have 
been  reported  in  speech  processing  literature  are  their  excellent  numerical 
behavior,  robustness  and  insensitivity  to  round-off  noise.  Ladder  forms  have 
also  been  used  very  successfully  for  adaptive  filtering  applications  like  noise 
cancelling,  and  channel  equalization. 

The  reasons  for  the  many  favorable  properties  of  ladder  forms  have  not 
been  fully  explored.  One  important  issue  is  the  interpretation  of  the  reflection 
coefficients  as  partial  correlations  between  the  forward  and  backward  innova¬ 
tions.  These  quantities  have  a  more  direct  physical  meaning  than  the  ARMA 
coefficients  {tij,  6*  J  .  which  are  related  in  a  complicated  way  to  the  poles  and 
zeroes  of  the  system  (through  finding  roots  of  a  polynomial).  Furthermore,  the 
mapping  from  ARMA  coefficients  to  poles  and  zeroes  is  very  sensitive,  i.e.,  small 
changes  in  the  coefficients  can  lead  to  relatively  large  changes  in  pole-zero  loca¬ 
tions.  The  sensitivity  to  changes  in  reflection  coefficients  seems  to  be  consider¬ 
ably  smaller.  Ladder  models  are  also  naturally  related  to  physical  models  like 
transmission  lines  or  a  layered  scattering  medium.  Such  models  correspond, 
for  example,  to  sound  waves  propagating  in  shallow  water.  The  reflection 
coefficients  correspond  to  actual  reflections  of  waves  propagating  in  the 
transmission  line.  The  {a*,  coefficients  do  not  seem  to  have  a  direct  physi¬ 
cal  interpretation  of  this  kind. 

In  view  of  these  observations  it  may  seem  somewhat  surprising  that  ladder 
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3.3  Event  Detection 

The  basic  concept  of  event  detection  and  its  usefulness  for  DSN  application 
was  discussed  already  in  [DSN].  Since  then,  we  performed  a  more  detailed  study 
of  the  ladder  form  as  an  event  detector,  in  the  context  of  speech  processing. 
The  following  is  a  summary  of  this  work,  as  described  in  [LM2]. 

The  many  nice  features  and  advantages,  such  as  a  fast  parameter  tracking 
and  excellent  convergence,  of  a  class  of  recursive  least-squares  ladder  estima¬ 
tion  algorithms  have  been  reported  by  us  in  [MLNV],  [MVL]  and  [ML2].  The  appli¬ 
cation  of  these  results  to  speech  modeling  and  synthesis  were  reported  in  [MLl]. 
In  [LM2]  we  present  a  novel  innovations  based  pitch  detection  technique  using 
the  ladder  algorithms.  The  important  features  of  this  pitch  detector  include: 

(i)  on-line  time-domain  operation  directly  on  the  speech  waveforms; 

(ii)  all  the  necessary  variables  for  the  pitch  detector  are  already  computed  by 
the  modeling  ladder  algorithms  and  therefore  is  well-suited  for  an 
integrated  implementation,  (e.g.,  VLSI). 

The  basic  assumption  is  that  the  speech  driving  process  consists  of  a 
approximately  Gaussian  part  (unvoiced),  and  a  jump  component  (voiced).  The 
pitch  positions  located  by  processing  the  innovations  alone  are  not  very  accu¬ 
rate  because  the  non-whiteness  of  the  Gaussian  component  will  act  to  cloud  the 
position  of  the  pulses,  mainly  due  to  phase-distortions.  In  our  ladder  algo¬ 
rithms,  a  gain  parameter  called  yp ,  is  computed  recursively  in  the  time-update 
recursions.  This  yp  parameter  turns  out  to  be  a  very  useful  Gaussian  log- 
likelihood  type  function  for  separating  the  mixed  driving  process  into  Gaussian 
and  non-Gaussian  components. 

Briefly,  the  yp  variable  appears  in  the  joint  Gaussian  distribution  of  speech 
samples  \yr-Vr-b  -Vr-pl  as 
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slopment  of  a  particular  class  of  filter  configurations,  the  ladder  or  lattice 
is,  which  allows  designs  to  meet  the  stringent  accuracy  and  speed  require- 
its  of  complex  applications  like  speech  synthesis  and  analysis.  Such  struc- 
»  go  beyond  the  realm  of  filtering  and  induce  novel  pipelined  and  parallel 
>rithms  for  the  solution  of  important  problems  in  matrix  arithmetic. 

In  this  paper,  we  first  describe  the  variety  of  operations  performed  by  an 
ptive  filter  for  speech  analysis.  A  short  review  of  the  CORDIC  algorithms 
irs  that  they  are  ideally  suited  to  realize  these  operations.  To  complete  the 
iy  of  the  filter  implementation,  a  real-time  speech  analysis  chip  is  discussed 
se  only  processing  elements  are  CORDIC  blocks. 

Next  we  turn  our  attention  to  more  general  numerical  problems  involving 
rices.  These  problems  are  classified  in  terms  of  the  order  of  the  matrices, 
eed  large  size  problems  generally  possess  some  structure,  for  instance 
rseness.  so  that  they  can  be  described  by  fewer  parameters  than  if  they  had 
arbitrary  number  of  elements.  Methods  for  the  solution  of  large  problems 
dictated  by  their  structure,  which  must  be  preserved  during  the  computa- 
ls  so  that  both  storage  and  execution  time  stay  low.  After  the  description  of 
bly  parallel  architectures  based  on  CORDIC  blocks  for  small  problems,  i.e., 
aout  any  structural  assumption,  we  survey  several  architectures  for  large 
blems.  Image  reconstruction,  beamforming  for  large  antenna  arrays  (with 
0  or  more  sensors),  seismic  data  processing,  boundary-value  problems  in 
tial  differential  equations  (such  as  the  ones  encountered  in  weather  forecast- 
1,  are  just  a  few  examples  of  large  problems  that  can  be  efficiently  solved  with 
se  architectures. 
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of  signal  processing  techniques  carries  over  to  numerical  linear  algebra:  direct 
techniques  are  related  to  triangular  matrix  factorization  and  inversion,  while 
transform  techniques  are  closer  to  orthogonalization  and  eigenvalue  determina¬ 
tion. 

We  address  the  issue  of  efficiently  computing  elementary  operations  which 
are  prominent  in  signal  processing  and  matrix  computations.  Most  signed  pro¬ 
cessing  structures,  whether  parallel  or  not,  tend  to  emphasize  the  need  for  fast 
multiplication  (e.g..  [4]).  However,  we  show  that  feist  vector  rotation  and  related 
coordinate  transformations  are  much  more  fundamental  elementary  operations 
for  these  problems.  Algorithms  that  perform  these  operations  were  introduced 
by  J.E.  Voider  [5],  who  named  them  CORDIC  for  Coordinate  Rotation  Digital  Com¬ 
puter.  Thus  we  advocate  the  use  of  CORDIC-like  blocks  as  elementary  proces¬ 
sors. 

Parallel  execution  is  needed  to  achieve  large  throughputs;  "highly  con¬ 
current"  systems,  consisting  of  many  processors  operating  simultaneously,  are 
developed  for  this  purpose.  Computer  scientists  have  designed  multiprocessor 
systems  and  array  processors  in  the  past,  see  e.g.  [3b]-[3d].  With  the  advent  of 
VLSI,  such  collections  of  processors  become  economically  feasible,  since  many 
processing  elements  may  be  realized  on  a  single  chip.  This  has  lead  to  much 
renewed  interest  in  parallel  processing  as  witness  [l]-[3].  Up  to  now,  the  effort 
has  been  mainly  to  restructure  for  parallel  execution  existing  algorithms, 
designed  for  execution  on  single  processor  machines.  However,  better  perfor¬ 
mance  can  be  expected  from  algorithms  designed  to  exploit  the  computational 
capabilities  of  parallel  processors. 

Much  can  be  learned  from  the  signal  processing  community,  which,  over  the 
last  two  decades,  has  paid  considerable  attention  to  issues  related  to  the  digital 
implementation  of  filters  or  transforms.  One  outcome  of  this  effort  has  been  the 
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resuited  in  a  push  for  the  development  of  faster  computing  structures  as  well  as 
algorithms  of  lower  computational  complexity.  For  particular  problems,  special 
purpose  hardware  has  also  been  developed,  but,  owing  to  its  lack  of  generality, 
this  approach  will  not  be  considered  here. 

General  purpose  uniprocessor  computers,  especially  microcomputers,  have 
been  utilized  in  the  high  speed  signal  processing  arena  with  only  limited  success 
for  primarily  three  reasons.  First,  they  cannot  in  general  compute  a  variety  of 
elementary  operations  such  as  multiplication,  vector  rotation  and  trigonometric 
functions  efficiently.  These  operations  Eire  very  common  in  signal  processing 
algorithms.  Secondly,  general  purpose  computer  architectures  provide  only 
cumbersome  address  arithmetic  for  data  structures,  such  as  circular  buffers, 
that  occur  frequently  in  communications  applications.  Finally,  signal  processing 
algorithms  exhibit  a  substantial  amount  of  parallelism  that  is  not  efficiently 
exploited  in  a  uniprocessor  system.  A  notable  exception  is  the  AMD2900  family 
which  allows  some  parallelism  through  the  extensive  use  of  two  port  random 
access  memories  (RAM’s). 

Signal  processing  techniques  roughly  fall  into  two  categories: 

direct  or  recursive  techniques,  where  the  signal  simply  passes  through  a 
digital  filter,  and 

-  transform  techniques,  where  the  signal  is  batch-processed.  A  global 
transformation  is  first  applied  to  a  batch  of  data,  then  simple  operations 
like  windowing  are  performed  in  the  transform  domain,  and  finally,  the 
resulting  signal  is  transformed  back  into  the  original  domain. 

Matrix  computations  share  many  features  with  digital  signal  processing. 
They  call  for  the  same  set  of  elementary  operations,  essentially  scaling,  rota¬ 
tions,  trigonometric  functions  and  square-roots.  Furthermore,  and  this  illus¬ 
trates  the  extent  of  the  connection  between  these  two  areas,  the  above  partition 
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more,  the  duality  between  Generalized  Levinson/Fast  Cholesky  algorithms  is 
clearly  displayed  as  a  construction/extraction  of  elementary  sections  of  the 
same  transfer  matrix  T, 

(hr)  A  VLSI  Speech  Analysis  Chip  Set  Based  on  the  Square-Root  Normalized 
Ladder  Form  [AULA] 

One  of  the  most  attractive  features  of  ladder  estimation  algorithms  is  their 
highly  modular  structure  and  suitability  for  VLSI  (Very  Large  Scale  Integration) 
implementations.  To  demonstrate  this  fact,  we  presented  in  [AMLA]  and  [AAM]  a 
chip  set  design  based  on  our  square-root  normalized  ladder  form.  It  is  shown 
that  the  equations  are  amenable  to  implementation  using  the  so-called  C0RD1C 
(Coordinate  Rotation  Digital  Computer)  algorithms  because  the  time  and  order 
updates  are  easily  representable  as  orthogonal  transformations  or  rotations.  An 
integrated  implementation  which  exploits  the  concurrency  of  the  ladder  recur¬ 
sions  together  with  possible  hardware,  speed  and  area  tradeoffs  are  presented. 
The  general  applicability  of  the  chip  set  to  other  signal  processing  tasks  is 
demonstrated  by  showing  that  the  discrete  Fourier  transform  (DFT)  is  naturally 
suited  to  the  architecture.  We  note  that  our  theoretical  development  of  the 
rotation  based  theory  for  ladder  forma  [MML]  [MMLD]  was  actually  motivated  by 
Implementation  issues  using  CORDIC'S. 

(v)  Highly  Concurrent  Computing  Structures  for  Digital  Signal  Processing  and 
Matrix  Arithmetic  [ADU] 

One  important  measure  of  the  utility  of  digital  signal  processing  algorithms 
has  traditionally  been  computational  complexity.  The  digital  implementation  of 
such  algorithms  has  been  frequently  difficult  (or  impossible)  since  they  tend  to 
be  compute  bound.  Consequently,  the  quest  for  real  time  signal  processing  has 
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(ii)  Hubert  Space  Array  Methods  for  Ladder  Realizations  [MML] 

In  this  paper  ire  presented  a  Hilbert  space  array  approach  for  deriving  fast 
estimation  algorithms  and  adaptive  signal  processing  algorithms,  that  are  recur¬ 
sive  in  time  and  order.  Ladder  (or  lattice)  forms  turn  out  to  be  the  natural  real¬ 
izations  of  these  algorithms.  From  a  stochastic  point  of  view,  the  natural  class 
of  processes  associated  with  these  techniques  include  stationary  but  also  non¬ 
stationary  processes.  They  sire  encountered  in  adaptive  signal  processing, 
speech  modeling  and  radar  and  sonar,  high-resolution  spectra]  estimation,  dis¬ 
tance  measures,  etc.  The  use  of  projections  and  orthonormalization,  e.g.,  via 
Gram-Schmidt  procedures  induces  real  and  complex  rotation  as  basis  opera¬ 
tions,  resulting  in  magnitude  normalized  variables  and  numericaUy  stable  com¬ 
putations. 


(iii)  £ -Contractive  Embeddings  of  the  Discrete  Lyapunov  Equation  [DGMVJ 

In  the  paper  [DMl],  a  Levinson-type  algorithm  for  fast  inversion  of  covari¬ 
ance  matrices  with  low  displacement  rank  has  been  derived.  In  [DGMV]  we  first 
present  a  normalized  version  of  these  recursions  (the  multichannel  case  is  con¬ 
sidered  throughout  the  paper).  Then,  the  quantities  involved  in  these  recursions 
are  embedded  into  a  state-space  description.  This  realization  exhibits  an  ortho¬ 
gonality  property  that  translates  into  the  para-unitarity  and  contractiveness  of 
its  associated  transfer  matrix  T(z). 

In  general,  a  para-unitary  transfer  matrix  can  be  shown  to  admit  a  cascade 
realization  with  sections  of  degree  one,  i.e.,  with  one  memory  element.  We 
demonstrate  that  the  Generalized  Levinson  algorithm  builds  recursively  a  cas¬ 
cade  realization  of  T(z),  where  the  elementary  sections  have  p  memory  ele¬ 
ments,  p  being  the  number  of  channels.  The  General  Fast  Cholesky  algorithm 
by  columns  [Ml]  is  also  shown  to  yield  the  same  cascade  realization.  Further- 
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certain  Riccatl-type  equations.  Indeed,  it  turns  out  that  the  state-space  system 
matrix  of  the  ladder  form  is  identical  to  a  certain  canonical  matrix  for  the 
discrete-time  Lyapunov  equation  called  Schwarz  form,  investigated  for  instance 
by  Mansour  [MAN],  in  the  context  of  Lyapunov  stability  test.  This  connection 
between  the  ladder  form  and  the  Schwarz  form  was  noted  by  Morf  [M2],  and  was 
elaborated  on  by  Lee  [LEE]. 

Since  Lyapunov  theory  is  intimately  related  to  well  known  classical  results 
of  the  polynomial  approach,  such  as  Routh-Hurwitz,  Hermite.  and  Schur-Cohn. 
etc.  the  ladder  form  solution  to  Lyapunov  and  Riccati-type  problems  often  offers 
much  simpler  and  elucidating  proofs  to  many  classical  results.  As  an  example, 
the  solution  to  the  classical  Schur-Cohn  problem  can  be  established  via  the 
ladder  form,  demonstrating  the  simplicity  of  the  ladder  canonical  form  method, 
compared  to  other  approaches  such  as  the  controller  form  solution  [BW]  and 
Hankel  matrix  approach  [DAT]. 

We  discuss  other  implications  of  the  state-space  structures  of  the  ladder 
form  such  as  demonstrating  the  extremely  efficient  way  of  generating  stationary 
states  of  any  white  noise  driven  system  in  the  minimal  number  of  steps  (equal  to 
the  order  of  the  system)  and  simultaneously  check  for  stability  or  stationarity  of 
a  process. 

We  may  note  that  representation  of  the  state-space  by  other  orthogonal 
polynomials  have  been  considered  in  the  literature.  For  example,  Good  [GO]  has 
considered  the  Chebyshev  polynomials  in  the  so-called  colleague  matrix  form. 
Barnett  [BAR],  and  Anderson  [AND]  have  considered  extensions  to  generalized 
orthogonal  polynomials  in  the  so-called  comrade  matrix  forms. 

We  briefly  describe  the  transformations  between  ladder  and  controller 
canonical  forms  and  outline  the  derivations  of  some  of  the  properties  mentioned 
above;  details  will  appear  in  a  forthcoming  paper. 
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Iheorotical  Studies 

In  addition  to  the  derivation  of  ladder  algorithms  of  different  types,  a  con* 
tiderable  amount  of  work  was  done  on  various  theoretical  issues  related  to 
ladder  structures.  As  mentioned  earlier,  ladder  forms  play  an  important  role  in 
system  theory.  In  fact  they  seem  to  provide  a  particularly  convenient  “coordi¬ 
nate  system"  for  studying  a  large  class  of  system  theoretic  problems.  We  briefly 
summarize  here  some  of  the  main  issues  that  were  addressed  so  far. 


(i)  State  Space  Structure  of  Ladder  Canonical  Forma  [HL3] 

The  purpose  of  this  paper  is  to  present  the  structure  of  the  ladder  canoni¬ 
cal  forms  in  a  state-space  model  context  and  to  recast  and  review  various 
results  involving  ladder  forms  from  a  state-space  model  view,  in  contrast  to  the 
input/output  (or  transfer  function)  point  of  view  stressed  in  the  past. 

The  transformations  between  ladder-  and  other  state-  space  canonical 
forms  are  presented.  By  means  of  state-space  models  we  established  some  of 
the  more  interesting  properties  of  the  ladder  forms,  such  as  being  the  natural 
canonical  form  for  the  discrete-time  Lyapunov  and  the  matrix  Riccati  equation 
for  stationary  and  so-called  nearly  stationary  or  finite  rank  processes. 

Conventional  canonical  state-space  realizations  are  obtained  by  construct¬ 
ing  bases  of  either  the  Hankel  matrix  of  the  Markov  parameters  or  of  the  con¬ 
trollability  (or  observability)  matrices  of  the  system.  The  ladder  form  realiza¬ 
tions  are  obtained  by  an  orthonormalization  (using  a  Gram-Schmidt  type  pro¬ 
cedure)  of  the  state-space  with  the  Szego  orthogonal  polynomials  as  basis 
[M1],[M2],  or  equivalently  via  (numerically  stable)  orthogonal  transformations. 
As  a  result  of  this  orthonormalization.  the  state  covariance  matrix  has  the  dis¬ 
tinct  advantage  of  being  diagonal  (or  unity  in  certain  normalized  case)  when  the 
input  is  white,  thus  making  it  an  attractive  choice  for  solution  to  Lyapunov  and 
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The  covariance  ladder  form  has  numerous  applications  in  adaptive  process¬ 
ing  and  it  also  provides  an  efficient  solution  technique  to  numerous  problems  in 
system  theory,  such  as: 

(1)  Minimal  realization  of  multi-input  and  multi-output  impulse  response 
sequences,  and  of  covariance  sequences  (the  so-called  stochastic  realiza¬ 
tion  problem); 

(ii)  Transforming  right  Matrix  Fraction  Descriptions  (MFD)  to  left  MFD's  and 
vice  versa.  This  problem  arises  in  the  ARMA  modeling  approach  to  the 
source  location  problem,  as  described  in  Section  2.4; 

(ill)  Model  approximation,  by  choosing  a  desired  number  of  sections  in  a  ladder 
realization; 

(iv)  Stability  tests  for  continuous  and  discrete  time-invariant  systems. 


Convergence  Analysis  [EM] 

An  important  step  in  the  development  of  any  new  recursive  algorithm  is  the 
analysis  of  its  convergence.  In  the  case  of  the  AR  ladder  form  such  an  analysis  is 
not  really  needed,  since  the  algorithm  provides  an  exact  solution  of  a  least- 
squares  problem,  and  the  asymptotic  properties  of  recursive  least-squares  esti¬ 
mators  are  by  now  well  known.  The  situation  is  very  different  in  the  case  of  the 
ARMA  ladder  which  involves  a  nonlinear  "bootstrapping"  technique.  In  [EM]  a 
preliminary  convergence  analysis  is  presented,  fhis  paper  gives  an  asymptotic 
analysis  of  a  particular  ladder  algorithm,  designed  for  auto-regressive  moving- 
average  (ARMA)  modeling  of  time-series.  A  thorough  analysis  of  a  second  order 
case  reveals  that  the  asymptotic  properties  can  be  improved  by  modifying  the 
algorithm.  The  modified  algorithm  is  then  shown  to  be  asymptotically  equivalent 
to  the  well-known  extended  least  squares  method  in  the  general  case. 
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This  would  be  the  case,  (or  example,  11  Uf  is  generated  by  linear  feedback  in  a 
control  system.  The  vector  process  zT  -  [yf,  uf]‘  will  then  have  an  AR 
representation 


p/r-i 

j^r-t 

*r-* 


(3.6) 


T 

By  applying  any  AR  modeling  technique  to  the  time-series  [zt  }J=i,  we  can  obtain 
estimates  of  the  ARMA  parameters  tq,  bt.  As  a  by-product  we  also  get  the 
parameters  c*,  dj  of  the  model  for  uj>.  If  these  are  not  needed,  it  is  possible  to 
simplify  the  algorithm  by  eliminating  the  equations  involving  /*,  p*. 


When  the  input  tty  is  not  known,  such  as  in  the  case  of  modeling  an  ARMA 
times  series,  the  situation  is  slightly  more  complicated.  The  unknown  input 
need  to  be  replaced  by  its  estimate,  which  under  certain  assumption  can  be 
shown  to  be  tr .  For  a  detailed  presentation  of  this  case  see  [FR2],  [LEE]. 

Similar  ideas  can  be  used  to  develop  ladder  forms  for  joint  process  estima¬ 
tion  [LHl],  [LEE]  and  for  ARMAX  models  [FR3]. 

The  square-root  normalized  ladder  algorithms  discussed  so  far  provide  a 
solution  to  the  so-called  pre-windowed  form  of  the  normal  equations.  Another 
version  of  the  normal  equations  that  is  commonly  used  in  speech  processing  and 
other  applications  is  the  covariance  form.  In  [PFM]  we  developed  the  growing 
memory  and  sliding  memory  covariance  ladder  algorithms,  using  a  simplified 
derivation  method.  Reference  [PFM]  also  includes  some  new  ladder  form  reali¬ 
zations  of  the  identified  models,  leading  to  convenient  methods  for  computing 
the  parameters  from  estimated  reflection  coefficients.  Another  topic  presented 
there  is  a  complete  solution  to  the  problem  of  possible  singularities  in  the 
ladder  update  equations. 
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greater  simplicity  of  the  normalized  recursions  leads  us  to  prefer  them  over 
their  unnormalized  counterparts.  Sometimes  it  is  desired  to  compute  the 
unnormalized  innovations  for  comparisons  with  other  spectral  estimation  tech¬ 
niques.  This  can  be  done  by  proper  scaling  of  the  normalized  variables: 

*p.r  -  (1  -  7p-i.r-iftRp^T  Sp.T  ~  unnormalized  innovations  (3.3a) 

Kps  =  (1  -  7P-i.r-i)*Ep,r  =  variance  normalized  innovations  (3.3b) 

where 

■« 

-  Rt  fl (1  -  Ki  j  K.t)*  ~  f  orward  innovations  covariance  (3.4a) 
(1  ~7j.-i.r-i)*  =  fl,  (1  -Kt-i  n.r-i)*  (3.4b) 

<*i 

Equations  (3.2a)-(3.2c)  represent  an  adaptive  whitening  filter.  Given  a  (vec¬ 
tor)  data  sequence  yp.  they  generate  an  innovations  sequence  tp  as  well  as  a 
set  of  reflection  coefficients  \Kpj\  which  fully  characterize  the  filter. 

Some  Extensions 

So  far  we  discussed  only  ladder  forms  for  multichannel  AR  processes.  As 
mentioned  before,  by  using  certain  embedding  techniques,  it  is  possible  to 
extend  these  results  to  more  general  situations.  In  [LEE],  [LFM]  and  [FR2]  we 
presented  such  an  extension  to  ARMA  models.  As  an  illustration  of  the  embed¬ 
ding  idea,  we  present  here  a  brief  description  of  its  application  to  the  ARMA 
case. 

Consider  the  problem  of  estimating  the  parameters  of  the  ARMA  model 

f 

given  by  Eq.  (3.1)  from  input/output  data  h/t.ut  {i*i-  Let  us  assume  that  the 
input  process  Uf  is  in  itself  an  ARMA  process: 

=  £  fi  Vr-i  -  £ 

«-l  <■! 


9i  uT^i  • 


(3.5) 


Thus,  with  just  one  extra  equation  we  can  solve  the  model  fitting  problem! 


Several  versions  of  the  exact  least  squares  ladder  forms  were  developed: 
the  unnormalized,  the  variance  normalized  and  the  magnitude  normalized 
ladder  forms.  The  normalizations  refer  to  the  signals  propagating  in  the  ladder 
filter,  i.e.,  the  forwards  and  backwards  innovations.  The  variance  normalized 
form  involves  signals  with  unit  variance  and  the  magnitude  normalized  form  has 
signals  with  magnitude  less  than  one.  The  simplest  of  these  three  forms  is  the 
magnitude  normalized  algorithm  which  is  summaried  below  for  the  all-pole  (AR) 
case.  It  should  be  noted  that  this  normalization  makes  it  possible  to  compute 
the  ladder  recursions  with  fixed-point  arithmetic,  which  is  ideally  suited  for 
micro-processor  and  very  large  scale  integration  (VLSI)  implementation! 

The  basic  square-root  normalized  AH  ladder  recursions  are  given  by 

Kp s  U  -  tp.T  tp.TfiKp+l.T-xU  -rv.T-irp.T-iY/z  +  tp.Trp.T- 1  (3.2a) 

Cp+I.r  *  [/  -Kp+u  f£+i'T]~*[tp.T  ~  Kp+i.Trp.T-i\  U  ~rp.T-\  tp.T-i)-*  (3.2b) 

rp*i.r  s  V  ~  Kp+i.r  ~  ^4-t.r  Bp.rl  (!  ~  Bj. t  Bp.r)"*  (3.2c) 

where 

Kf  =  reflection  coefficient  matrix 

tp  -  magnitude  normalized  forward  innovations 

rp  s  magnitude  normalized  backward  innovations 


The  notation  [  •  f*  denotes  the  matrix  square  root.  The  derivation  of  the  all- 
pole  normalized  ladder  forms  can  be  found  in  [LEG],  [LMl].  Alternatively,  one 
may  use  the  unnormalized  ladder  forms  [MLl]  [VT.r>l.  However,  preliminary 
tests  seem  to  indicate  that  the  normalized  recursions  perform  at  least  as  well  as 
(and  maybe  better  than)  the  unnormalized  version.  This  fact  combined  with  the 
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windowing  (e.g.,  as  in  the  covariance. or  auto-correlation  method).  An  assump¬ 
tion  of  stationarity  (Toeplitz  covariance  matrix)  is  used  in  these  computations. 

Block  processing  is  often  useful,  but  in  adaptive  processing  applications 
where  computations  have  to  be  carried  out  in  real-time,  recursive  methods  are 
needed.  Recursive  methods  have  distinct  advantages  even  in  situations  where 
block  processing  could  be  used.  They  are  generally  easier  to  implement  and 
lead  to  greater  flexibility.  In  signal  processing  applications  where  the  signals 
have  time-varying  nonstationary  statistics  it  is  useful  to  be  able  to  update  con¬ 
tinuously  the  estimates  of  the  model  parameters  (e.g.,  the  reflection  coefficients 
of  the  ladder  form),  rather  than  do  so  at  fixed  time  intervals.  This  allows  faster 
adaptation  of  the  filter  to  variations  in  the  data.  Because  of  these  and  other  rea¬ 
sons  there  is  great  interest  in  developing  recursive  algorithms  for  fitting  ladder 
filters  to  an  observed  time  series.  Some  results  have  been  reported,  using  gra¬ 
dient  search  techniques  which  adaptively  vary  the  reflection  coefficients  so  as  to 
minimize  a  square  error  criterion.  These  gradient  techniques  often  show  slow 
convergence  which  limits  their  usefulness. 

The  exact  least-squares  method  is  a  recursive  technique  for  computing  the 
ladder  coefficients.  It  computes  at  each  time  step,  a  set  of  reflection 
coefficients  which  minimizes  the  sum  of  the  squared  errors  for  the  data 
sequences  up  that  time.  Thus,  it  provides  an  exact  solution  to  the  least  squares 
model  fitting  problem!  The  computation  is  recursive  both  in  time  and  in  model 
order,  which  adds  an  extra  degree  of  flexibility.  What  is  rather  remarkable  is 
that  this  exact  recursive  solution  is  obtained  at  almost  no  increase  in  the  com¬ 
putational  requirements.  The  normalized  ladder  forms  involve  only  three 
update  equations  per  time  step  and  per  section  of  the  filter.  Two  of  these  equa¬ 
tions  are  needed  to  implement  the  ladder  filter  (i.e.,  propagate  the  forwards  the 
backwards  innovations),  while  the  third  updates  the  reflection  coefficients. 


-36- 


forms  ve  not  more  widely  used  in  adaptive  processing.  One  reason  for  this  may 
be  the  fact  that  until  fairly  recently  techniques  for  fitting  ladder  structures  to 
observed  data  were  restricted  to  very  simple  classes  of  systems  of  the  single¬ 
input  single-output  all-pole  (or  all-zero)  type.  Many  adaptive  processing  prob¬ 
lems  involve  ARMA  models  and  multiple  inputs  and  outputs.  Using  a  certain 
embedding  approach  we  were  able  to  extend  considerably  the  class  of  problems 
to  which  ladder  structures  can  be  applied. 

A  second  reason  which  may  have  limited  the  use  or  ladder  structures  was 
the  lack  of  efficient  recursive  computational  methods  for  estimating  ladder 
parameters  (reflection  coefficients).  Until  fairly  recently,  the  only  techniques 
available  for  ladder  modeling  were  of  the  block  processing  type,  i.e.,  they 
operated  on  data  collected  over  a  finite  time  interval.  Adaptive  processing 
requires  frequent  updating  of  the  estimated  parameters  which  is  best  achieved 
by  a  recursive  algorithm.  It  is  possible  to  use  block  processing  algorithms  in  a 
sliding-window  mode  to  obtain  point-by-point  updates,  but  this  often  requires 
prohibitive  amounts  of  computation.  The  development  of  the  exact  least 
squares  recursions  for  ladder  forms  provides  a  strikingly  simple  and  efficient 
algorithm,  which  makes  the  use  of  ladder  forms  in  adaptive  processing  much 
more  attractive. 

Recursive  Ladder  Farms 

Various  algorithms  have  been  proposed  for  finding  a  ladder  form  represen¬ 
tation  of  time  series.  These  algorithms  are  mostly  of  the  block-processing  type. 
Typically,  the  algorithm  is  based  on  the  fact  that  given  the  data  covariance 
matrix,  the  reflection  coefficients  can  be  computed  by  a  Levinson-type  recur¬ 
sion.  Since  the  data  covariance  is  usually  unknown,  it  has  to  be  replaced  by  an 
estimate.  This  is  done  by  using  the  sample  covariance  of  the  data  with  a  proper 
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P(Vt . VT-p)  =  l2TrRp|“1/8exp{- |-7p.rl  .  (3-7) 

where 

7p.T  -  [ VT . Vr-p]  Rp.V  [ Vt . VT-pY  (3-8) 

and  R p  is  the  covariance  matrix.  Taking  the  logarithm,  we  have  a  log-likelihood 
function  of  the  form 

Ap  =  In  |RP  |  +  llyll*., 

=  In  |Rp  |  +  .  (3.9) 


It  turns  out  that  the  increase  in  Ap  per  sample,  i.e.,  the  derivative,  6  Ap  of 
Ap  ,  is  a  very  sensitive  measure  of  the  "unexpectedness”  of  a  speech  sample, 
i.e.,  a  measure  of  the  deviation  of  the  actual  distribution  from  the  Gaussian 
hypothesis.  This  is  actually  a  very  useful  approach  in  the  so-called  innovations 
based  detection  of  outliers  in  the  statistics  literature,  see  [ANS],  [AN]. 

Thus  a  very  simple  but  powerful  pitch  pulse  detector  can  be  designed  as  fol¬ 
lows: 

(i)  use  a  local  maximum  algorithm  either  on  6  Ap,  or  Syp,  returning  a  zero 
value  if  the  maximum  is  below  some  threshold; 

(ii)  multiply  the  resulting  log-likelihood  by  the  innovations,  tp ,  and  then  use  a 
simple  exponential  threshold  detector  similar  to  that  of  [GR]  to  locate  the 
pitch  pulses. 

A  sample  of  the  pitch  detection  results  is  shown  in  the  sequence  of  figures 
below. 
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Figure  a4.  10th  Order  Innovations 
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Figure  3.6.  6  yio  (Log-likelihood  Ratio) 
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figure  3.8.  Detected  Pitch  Pulses 
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We  may  note  that  this  particular  pitch  detection  scheme  was  designed  to 
illustrate  the  basic  principles  and  the  power  of  our  ladder  form  algorithms,  and 
is  by  no  means  an  optimized  solution  to  the  problem  of  pitch  detection.  This  is 
clear  from  the  fact  that  plosive  sounds  also  generate  large  7's;  however  they 
tend  to  be  isolated  and  not  to  occur  in  periodic  intervals. 

From  a  stochastic  modeling  point  of  view,  one  should  be  able  to  decompose 
the  driving  process  into  essentially  continuous  and  discontinuous  or  jump  type 
components.  We  view  the  above  scheme  as  a  first  step  in  this  direction. 

The  algorithm  discussed  above  can  be  used  as  a  sensitive  detector  of  vari¬ 
ous  events  in  a  DSN  system  such  as:  the  appearance  (or  disappearance)  of  a  new 
source  in  the  area  under  surveillance;  a  change  of  course,  in  the  flight  path  of  an 
airplane  or  other  source  of  acoustic  energy;  a  change  in  the  spectral  charac¬ 
teristics  of  the  source.  This  event  detector  can  be  used  to  identify  which  por¬ 
tions  of  the  large  amount  of  data  that  are  expected  to  be  available  in  a  DSN  sys¬ 
tem,  are  of  interest.  We  then  can  proceed  to  process  in  detail  only  the  selected 
data,  thus  achieving  considerable  savings  in  computational  requirements. 

We  note  that  other  types  popular  log-likelihood  functions  distance  or  distor¬ 
tion  measures  can  be  computed  efficiently  via  ladder  forms,  for  a  discussion  of 
such  results  see  [LM4], 

Application  to  High  Speed  Modem  and  Packet  Radar/Sonar 

Distributed  Sensor  Networks  clearly  involve  not  only  collecting  and  process¬ 
ings  data,  but  communication/exchange  of  digital  information  as  well.  We  also 
stated  earlier,  that  event  detection  was  an  Important  component  of  data  pro¬ 
cessing.  Events,  by  definition,  are  discrete  or  discontinuous  phenomena  that 
generally  require  different  sets  of  mathematical  tools  and  processing  algorithms 
than  continuous  or  "slowly''  varying  objects,  such  as  spectral  parameters  or 
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track  information  of  slowly  moving  sources.  Digital  information,  is  discrete  by 
definition,  hence  a  digital  communication  system  has  to  be  based  on  signals  that 
exhibit  discontinuous  components  or  events. 

Our  work  on  tools  for  event  detection,  e.g.,  using  ladder  forms  for  comput¬ 
ing  likelihood  functions,  is  now  clearly  applicable  to  processing  digital  communi¬ 
cation  data.  Although  such  work  has  already  been  done  in  the  area  of  digital 
communication  in  general  and  modems  (modulators  -  demodulators)  in  particu¬ 
lar,  relatively  few  results  are  available  on  the  computational  complexity  and 
implementation  of  optimal  modem  (e.g..  9.6  kbits  for  dial-up  lines)  designs.  We 
only  note  that  currently  available  modems  are  either  far  from  the  telephone 
channel  capacity  (e.g.,  30  kbits  for  a  single  commercial  dial-up  telephone  chan¬ 
nel)  or  very  expensive.  The  availability  of  high  capacity  digital  links  is.  however, 
a  crucial  parameter  in  a  high  performance  DSN. 

Our  work  on  fast  computation  of  likelihood  functions  and  other  measures 
using  ladder  forms  has  now  been  extended  to  include  the  problem  of  optimal 
design  of  digital  modems.  As  an  example  that  would  test  our  mathematical 
design  tools,  we  decided  to  work  out  a  design  of  binary  FSK  and  PSK  moderns. 
This  design  is  based  on  ladderforms,  hence  an  economic  implementation  could 
be  based  on  our  LADMOS  -  VLSI  chip  design. 

Our  preliminary  simulation  studies  have  indicated  that  this  modem  has 
Bayes  optimal  detection  characteristics  in  colored  Gaussian  noise. 

In  a  DSN  context,  the  combination  of  passive  or  active  sensors  (e.g., 
sonar  /radar  /acoustic),  and  analog  or  digital  communication  is  an  extremely 
promising  solution,  in  our  opinion,  to  the  design  constraints  imposed  on  a  DSN 
by  its  distributed  character  and  most  likely  limited  communication  capacity 
over  hard  channels  (cables,  fiber  optic,  etc.). 
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We  envision  that  cooperative  sensors  and  possibly  active  transmitters  for 
both  source  tracking  and  communication  will  lead  to  new  types  of  solutions  for 
DSN  problems.  For  example,  the  bistatic  radar  mode  of  operation  is  very 
natural  between  communicating  sensor  nodes,  if  transmitters  are  linked  to  the 
sensors;  hence  cooperative  tracking  of  sources,  as  well  as  "space-hopping”  (in 
addition  to  frequency  and  time  hopping  radar  /sonar  modes)  is  clearly  achiev¬ 
able.  In  addition  if  transmitting  signals  are  modulated  with  digital  data,  high 
bandwidth  digital  communication  channels  become  available  between  sensor 
nodes!  I.e.,  we  can  envision  packet  radar  /sonar  networks.  In  order  to  achieve 
this  goal  we  clearly  need  combined  event  detection  and  digital  modem  related 
signal  processing  tools  of  the  type  we  are  developing. 
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3.4  Doubling  and  Tree  Algorithms  [M3],  [MDl],  [ADM] 

The  inversion  of  Toeplitz  covariance  matrices  arises  in  practically  all  esti¬ 
mation  problems  related  to  stationary  time  series.  If  \yt  j  is  a  stationary  ran¬ 
dom  process,  and  *  is  some  statistically  related  random  variable,  then  the 
least  squares  estimate  of  x  based  on  observations  \yt  i«*o  is  given  by 

x  =  .  (3.10) 

where  Ry  is  the  covariance  matrix  of  the  data  vector  Y  -  [y’0,  ....  y'r\  and 
Rya  is  the  cross  covariance,  i.e.,  Ry,  ~  E  {yx\.  Because  of  the  assumption  of 
stationarity,  Ry  has  a  Toeplitz  structure.  The  random  variables  [yt  {/Lo  do  not 
necessarily  represent  a  time  sequence.  They  may  be,  for  example,  measure¬ 
ments  at  the  outputs  of  multiple  sensors,  in  which  case  the  index  t  represents 
the  senses'  number  (i.e.,  it  is  a  spatial  rather  than  a  temporal  index).  This  type 
of  problem  arises  in  adaptive  beamforming,  adaptive  antenna  arrays,  and  in  a 
DSN  system  consisting  of  multiple  acoustic  or  electromagnetic  sensors. 
Development  of  efficient  solution  techniques  for  Toeplitz  sets  of  equations  is 
therefore  of  great  interest  in  the  context  of  DSN  signal  processing. 

Using  the  HGCD  algorithm  of  Abo.  Hopcroft,  UUman  [AHU],  Gustavson  and 
Yun  [GY]  constructed  an  algorithm  to  invert  a  Toeplitz  matrix  in  0(n  log®n). 
The  main  principle  used  in  the  HGCD  (half -greatest-common-divisor)  as  well  as  in 
our  new  algorithm  Is  the  "divide  and  conquer”  or  "doubling"  approach.  In  both 
cases  an  implicit  form  of  the  inverse  of  an  n  x  n  Toeplitz  matrix  T  is  found  in 
0(n  log^)  operations;  however,  our  algorithms  can  work  with  a  larger  class  of 
matrices,  and  they  have  the  potential  for  being  more  efficient,  e.g.,  0{n  log  n). 

We  are  interested  in  using  the  idea  of  doubling  directly  on  Toeplitz  matrices 
without  referring  to  the  HGCD  algorithm  and  to  find  expressions  for  inverting  a 
more  general  class  of  matrices,  that  is  related  to  Toeplitz  matrices.  The  objec- 
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tive  was  to  reduce  the  0(n®)  operations  of  the  well-known  algorithm  of  Levinson 
[LEV]  and  its  a-Toeplitz  matrix  extensions  [Ml],  [FMKL],  [KKM]. 

First  we  recall  the  definition  and  properties  of  displacement  ranks  required 
in  the  proof  of  the  basic  algorithm.  Then  we  introduce  the  Fast  Toeplitz  Matrix 
Inversion  Algorithm.  It  is  based  on  partitioning  of  a  matrix 

=  [£*]•  <311> 

where  the  indices  indicate  the  dimensions  of  matrices,  and  its  inverse 

rv  rw 

*2  =  Cn  f  n  '  (312) 

The  Schur  complements  of  the  submatrices  At  and  Dn  are  defined  as 

An  -  An  -  Bn  AT1  Cn  ,  A.  =  A.  -  A.  AT1  A,  .  (3.13) 

Now  we  have  the  following  well  known  relations 


Aj»  =  An 1  s  Art*  +  Aw*  At  A  n  At  At 1 

(3.14) 

Bn  5  -An  1  At  B n  —  ~An  Bn  Dn  1 

(3.15) 

C n  ®  “"An  C,,  At 1  =  -At 1  Cn  An 

(3.16) 

=  Bn1  =  AT1  +  AT1  At  An  At  AT1 

(3.17) 

The  algorithm  HI!  (half  a-Toeplitz  matrix  inversion)  follows  by  applying  the 
properties  of  Toeplitz  and  sums  of  products  of  Toeplitz  matrices  to  these  formu¬ 
las,  and  calling  the  HTI  algorithm  recursively  for  the  necessary  inversions  of 
submatrices.  The  crucial  observation  is  that  the  Schur  complements  have  an 
efficient  representation  that  is  just  as  complex  as  the  submatrices  themselves, 
Le..  they  have  the  same  so-called  displacement  rank,  a  number  that  measures 
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for  multiplications  and  inversions. 

We  conclude  by  mentioning  some  applications  of  these  results.  Matrix 
inversion  algorithms  that  are  based  on  partitioning  lead  to  the  problem  of 
inverting  the  Schur  complement  of  a  submatrix.  If  one  is  interested  in  algo¬ 
rithms  that  can  take  advantage  of  the  structure  of  the  matrix  to  be  inverted, 
the  question  arises,  what  properties  are  invariant  under  the  Schur  complement 
action.  As  an  illustration,  we  sketch  a  partial  list  of  matrices  that  can  be  shown 
to  have  such  invariance; 

Hankel  and  Hankel  plus  Toeplitz  type  equations; 

Banded  and  rational  (ratio  of  banded)  type  matrices,  0(n  log  6); 

Block  Toeplitz  circulant  or  Hankel  matrices; 

Matrices  that  are  related  to  two-dimensional  (2-D)  and  multi-dimensional 

(M-D)  problems,  e.g.,  Toeplitz  block  Toeplitz  matrices; 

Other  combinations  of  the  above. 

There  are  various  other  interesting  topics  that  are  related  to  doubling  algo¬ 
rithms  of  this  type,  e.g.. 

Stochastic  Interpretation  of  doubling /halving  algorithms  [DMl],  [MDl]; 

Fractal  ladder  form  realization  of  matrix  inversion. 

open  problems,  especially  computational  refinements, 

other  applications,  estimation,  detection,  pattern  recognition, 

-  image  reconstruction  from  projection,  etc. 

The  study  of  structures  of  sets  that  are  inherited  in  subset  is  a  topic  by 
itself.  Generally  this  leads  to  recursive  function  theory.  A  very  interesting  class 
of  objects,  perhaps  up  to  now  more  of  a  curiosity,  consists  of  the  so-called- frac¬ 
tals.  They  are  best  thought  of  as  objects  created  by  recursive  function  calls  or 
by  "mirror  galleries".  For  an  example  see  the  "fractal  ladder  realization"  in 
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4] 

Many,  especially  yet  undiscovered  applications  can  be  obtained  by  suitable 
orations  that  convert  structure  of  matrices  into  sparseness.  Simple  examples 
this  type  are  matrices  that  have  higher  order  displacement  rank,  e.g.,  [KKM]. 
ich  concepts  are  useful  for  matrices  with  n 01 -order  polynomial  variation  along 
Btgonals  and  off-diagonals,  compared  to  constants  for  Toeplitz  matrices. 
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3.5  Adaptive  Signal  Proc easing 

The  distributed  sensor  network  operates  in  an  uncertain  environment. 
Various  system  parameters  are  unknown  to  the  processing  nodes  either  because 
of  random  effects  or  because  of  their  time  varying  nature.  Therefore,  fixed  sig¬ 
nal  processing  schemes  may  be  inadequate  in  certain  DSN  applications,  and  an 
adaptive  processing  technique  will  have  to  be  used.  Several  adaptive  techniques 
were  developed  as  part  of  our  research  program,  and  are  briefly  summarized 
below. 

Adaptive  linearfluse  Filters  [FH2] 

An  important  class  of  filters  commonly  used  in  digital  signal  processing  is 
the  class  of  finite  impulse  response  (FIR)  linear-phase  filters.  Such  filters  are 
needed  in  applications  where  frequency  dispersion  due  to  nonlinear  phase  is 
harmful,  e.g.,  speech  processing  and  data  transmission.  Signals  in  the  passband 
of  linear-phase  filters  are  reproduced  exactly  at  the  filter  output,  except  for  a 
delay  corresponding  to  the  slope  of  the  phase  vs.  frequency  plot.  An  extensive 
literature  exists  on  the  design  of  linear-phase  filters  with  a  desired  frequency 
response  and  various  efficient  design  algorithms  have  been  developed.  Most  of 
this  literature  deals,  however,  with  non-adaptive  processing,  i.e.,  designing  the 
filter  to  have  a  known  frequency  response  (or  Impulse  response).  Relatively  lit¬ 
tle  work  seems  to  have  been  done  on  adaptive  linear- phase  filters.  i.e.,  filters 
whose  characteristics  are  determined  on  the  basis  of  an  observed  time-series, 
and  not  on  apriori  specifications. 

An  FIR  linearphase  filter  is  characterized  by  a  symmetric  impulse 
response.  Let  A(s)  represent  the  transfer  function  of  am  FIR  filter,  where 

A(s)  =  o0  +  a,*-1  +  •••  +  ans“n  (3.18) 
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z~l  =  unit  delay  operation,  i.e. ,  z~lxt  =  x4_i  . 


his  filter  will  be  a  linear-phase  filter  if  a*  =  a*-*.  In  this  paper  we  will 
s  assume  that  a0  =  1.  Such  filters  appear  in  adaptive  processing  applica- 
often  without  being  explicity  identified  as  having  linear  phase.  As  an 
pie  consider  the  high  resolution  spectral  estimation  of  narrowband  signals, 
toregressive  (AR)  modeling.  Various  spectral  estimation  techniques  such 
s  Maximum  Entropy  method  [BUR],  or  Pisarenko's  method  [PIS],  involve 
:ting  of  as  AR  predictor  A(z)  to  the  observed  data  sequence  yt  and  then 
it  to  compute  the  spectrum 


Sy(o)  - 


A(z)A(z~') 


!*=•>« 


(3.19) 


’  the  data  yt  has  a  pure  line  spectrum  (i.e.,  a  sum  of  sinusoids  of  different 
sncies)  it  is  straightforward  to  show  that  A(z)  will  be  a  symmetric  polyno- 
with  all  its  roots  on  the  unit  circle.  In  other  words,  the  predictor  A(z)  is 
cial  case  of  a  linear-phase  filter.*  In  many  of  these  techniques  the  sym- 
of  A(z)  is  imposed  implicitly. 

be  interconnection  between  the  predictor  of  a  narrowband  signal  and 
-phase  filtering  should  make  the  techniques  described  in  this  paper  useful 
plications  including  adaptive  line  enhancement,  adaptive  noise  cancelling 
laptive  array  processing.  Another  area  where  linear-phase  filters  are  of 
1  interest  is  adaptive  channel  equalization. 


le  objective  of  this  paper  is  to  formulate  the  adaptive  linear-phase  fllter- 
oblem.  discuss  its  interpretations  in  the  context  of  linear  least-squares 
tion  and  derive  several  types  of  adaptive  algorithms.  The  emphasis  is  on 
oeral  features  of  the  adaptive  linear-phase  filter  and  not  on  any  specific 

«ral,  *«  is  a  zero  of  a  Hnear^phaae  filter,  so  is  l/z(.  For  a  zero  on  Ih  unit.dr- 
—  •  *  ,  and  therefore  1/*,  =  S-’**®  .  The  complex  pair  je,U®,e 
ts  a  sinusoid  of  frequency  Ua. 
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applications.  While  many  aspects  of  these  adaptive  filters  have  been  discussed 
in  the  literature,  no  unified  treatment  of  this  subject  seems  to  be  available. 

Adaptive  line  Enhancement  [NH] 

Recently  there  hats  been  great  interest  in  the  problem  of  separating  narrow 
band  and  periodic  signals  from  wide  band  signals  or  noise,  using  least-squares 
predictors.  A  considerable  amount  of  work  was  based  on  gradient  type  tech¬ 
niques,  referred  to  as  Adaptive  Line  Enhancer  (ALE),  by  B.  Widrow  et  al.  [WG].  In 
[NM]  we  analyze  the  optimal  behavior  of  the  ALE,  which  coincides  with  the  exact 
least-squares  predictor  in  steady  state.  The  ALE  and  the  noise  canceller  can  be 
viewed  as  a  special  case  of  the  joint  estimator  process  estimator  or  instrumental 
variable  method,  see  e.g..  [Ml 2]  and  [FR4],  where  the  reference  input  is  a 
delayed  (or  suitably  pre-flltered)  version  of  the  primary  input,  as  shown  in  Fig¬ 
ure  3.9. 


Figure  3.9.  Sock  Diagram  of  the  ALE 
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A  more  detailed  examination  led  us  to  define  several  different  but  related 
itions  of  displacement  rank,  leading  to  a  classification  of  algorithms  and  to  the 
testion  of  minimality  of  the  various  associated  representations. 

In  the  discrete  case,  a  new  Toeplitz  block  matrix  embedding  principle  pro- 
des  us  not  only  with  an  elegant  framework  for  relating  various  algorithms  for 
living  Toeplitz  systems,  but  also  with  an  interpretation  of  the  embedding 
atrix  as  a  joint  covariance  matrix.  Because  of  the  larger  size  of  this  matrix, 
ie  number  of  random  variables  is  roughly  alpha  times  the  original  number  of 
iriables,  hence  we  encounter  an  alpha  multiplicity  in  a  process  representation. 

This  leads  to  the  first  physical  Insight  that  we  can  now  provide:  alpha  is  a 
lultiplicity  in  the  representation  of  a  process,  as  the  sum  of  the  outputs  of 
Ipha  linear  time-invariant  (but  not  necessarily  finite-dimensional)  systems 
Larting  at  time  zero,  i.e.  the  superposition  of  their  transients.  An  alternative 
iterpretation  turned  out  to  be  to  model  the  process  as  the  prediction  error  of 
ne  process  given  alpha  other  ’'pseudo''  processes;  pseudo  processes  because 
he  estimation  errors  are  by  definition  orthogonal  to  the  "given"  processes,  but 
his  implies  that  these  "given"  processes  are  not  observable  in  the  process  to  be 
lodeled  by  the  projection  theorem.  The  alpha  underlying  innovations  in  the 
given"  processes  can  now  be  viewed  as  "pseudo”  innovations.  Like  the  innova- 
Lons  they  are  white  processes  that  can  generate  the  process  causally;  however, 
dike  the  true  innovations  they  can  in  general  not  be  obtained  In  a  causal  way, 
e.  they  are  smoothing  residuals,  and  they  have  in  general  an  unobservable 
omponent.  This  component  can  be  viewed  as  an  additive  "dither"  noise  or 
enciphering  key”  that  allows  the  alpha  pseudo  innovations  to  be  white.  In  par- 
Icular  the  last  reference  points  out  some  potential  applications  of  this 
epresentation. 

A  careful  analysis  shows  that  in  general  a  mixed  representation  can  be 
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sumption  of  stationarity  is  very  widely  used.  One  reason  for  this  dicho- 
ippear  to  be  observations  of  abrupt  transition  between  the  complexity  of 
lary  and  non-stationary  results.  We  have  recently  introduced  an  index  of 
tarity  (alpha),  [FMKL]  that  has  the  potential  of  providing  a  more  gradual 
action  of  non-  stationarities  into  the  modeling  of  stochastic  processes. 

le  significance  of  this  index  alpha  is  that  the  larger  the  value  of  alpha  the 
ion-stationary  a  process  is.  In  addition,  the  complexity  of  models  and  cal- 
>ns  increases  only  proportional  to  this  index,  in  a  similar  way  as  the  order 
ndex  of  complexity  for  finite  dimensional  linear  systems.  However,  there 
■ocesses  with  finite  alpha  and  infinite  order  (i.e.,  non-rational)  and  vice 
hence  the  two  notions  are  quite  independent,  except  for  constant  param- 
taite  dimensioned  linear  systems  where  alpha  is  bounded  above  by  the 
[KA2],  [MSK].  The  difference  is  seen  from  the  fact  that  finite  order  is 
lent  with  semi-separability  of  the  covariance  kernel,  say  whereas 

is  given  by  the  rank  of  the  separable  displacement  function 
,<)  =  dK{s,t)/dt  +  dK(s,t)/  ds  .  see  e.g.,  [KLM],  where  _|  is  an 
or  of  the  divergence  type.  Livshits  [LY]  defined  an  equivalent  infinitesimal 
ition  function,  and  its  separability  rank  the  rank  of  the  process.  For 
te  processes  we  defined  in  [Ml]  the  "shift"  or  tensor  rank  for  matrices 
jo  Toeplitz  matrices,  or  what  we  call  now  the  displacement  rank  -  alpha, 
matrices  are  basically  rational  functions  of  Toeplitz  matrices,  whereas 
>rder  systems  lead  to  matrices  of  rational  functions. 

my  estimation,  communication  and  control  problems  require  the  solution 
plitz  and  related  systems  of  equations  or  integral  equations  of  the  dis¬ 
tent  type.  Such  solutions  can  be  obtained  efficiently  via  the  use  of  exist- 
tthods,  such  as  the  Krein-Levinson  type  algorithms  [KLM],  and  the  more 
doubling  type  algorithms  [M3]. 
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nels  of  Livsic,  [L]. 

The  displacement  rank  used  so  far  is  actually  called  a  "strict"  displacement 
rank  in  contrast  with  two  other  displacement  ranks,  the  so-called  "+"  and 
displacement  ranks,  see  e.g.  [KKM].  All  these  indices  are  closely  related;  their 
difference  is  bounded  by  two.  It  turns  out  that  the  ChristofTel-Darboux  type 
representation  of  the  inverse  obtained  using  the  strict  displacement  rank 
[FMKL]  is  not  always  minimal.  However  it  is  possible,  applying  the  semi-group 
property  of  Schur  complementation,  see  [A],  to  find  a  minimal  representation  of 
the  inverse  that  can  be  exploited  to  obtain  still  another  Levinson-type  algorithm 
far  inverting  matrices  close  to  Toeplitz.  The  interest  of  this  new  algorithm  lies 
more  in  the  insight  it  gives  on  stochastic  process  representation  and  ladder 
forms  than  into  the  computational  savings  which  are  minor.  However  the  same 
type  of  minimal  representation  can  also  be  obtained  for  quantities  of  interest 
(e.g.  Schur  complement)  in  the  doubling  algorithm  described  in  [M3]  and  the 
computational  complexity  is  significantly  reduced  there.  The  quantities  involved 
in  the  minimal  Levinson-type  algorithm  also  appear  very  closely  related  to  the 
ones  in  the  doubling  algorithm. 

The  notion  of  displacement  rank  and  the  associated  Toeplitz  algebra  are 
known  to  extend  to  operators,  i.e.,  continuous  time  models,  see  e.g.  [KLM]. 
Furthermore  the  Generalized  Levinson  equations  for  nondisplacement  kernels 
were  derived  in  that  paper.  In  this  last  section  a  continuous  time  version  of  the 
doubling  method  [M3]  will  be  described  that  yields  yet  another  efficient  compu¬ 
tational  procedure  for  determining  Fredholm  resolvents  of  nondisplacement 
kernels. 

Ftaudo  Innovations  Representations  for  Alpha-Stationary  Processes  [M DB] 

Almost  all  processes  found  in  practical  situations  are  non-stationary,  yet 
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cks  of  different  or  varying  sizes,  a  structure  of  special  interest  in  state-space 
dels  and  minimal  realizations. 

Various  embeddings  of  matrices  close  to  Toeplitz  in  the  sense  of  [Ml],  Le., 
lple  rational  functions  of  Toeplitz  matrices,  can  be  used.  Two  useful  Toeplitz 
ck  embeddings  for  symmetric  matrices  are  inferred  from  their  representa- 
n  either  as  a  Toeplitz  matrix  plus  an  algebraic  sum  of  strictly  lower  triangular 
split z  matrices  times  their  transpose  or  a  Toeplitz  matrix  minus  a  sum  of 
rer  and  upper  triangular  Toeplitz  matrices  times  their  transpose  (mixed 
>resentation).  Application  of  the  Toeplitz  block  version  of  the  LWR  algorithm 
the  first  embedding  matrix  provides  a  simple  derivation  of  the  Generalized 
rtnson  recursions  of  [Ml],  [FMKL]  while  the  second  embedding  yields  a 
rtnson-type  algorithm  with  improved  numerical  properties  since  the  embed- 
ig  matrix  is  positive  definite  when  the  original  matrix  is  positive  definite. 

Consider  the  case  where  the  original  matrix  is  the  covariance  matrix  of  a 
icrete-time  stochastic  process.  The  two  Toeplitz  block  matrix  embeddings 
esented  above  have  (r+l)x(r+l)  blocks  where  r  is  the  displacement  or  shift 
ok  of  the  covariance  matrix,  an  index  of  non-stationarity  of  the  process  [Ml], 
ire  over  the  embedding  matrix  associated  with  the  mixed  representation  can 
interpreted  as  a  joint  covariance  matrix  A  non-stationary  process  may 
erefore  be  modeled  as  the  prediction  error  of  one  process  given  r  other 
■eudo"  processes.  The  r  underlying  innovation  processes  are  indeed  "pseudo" 
novations;  like  the  innovations  they  are  white  processes  that  can  generate  the 
ocess  causally,  however,  unlike  the  true  innovations  they  cannot  in  general  be 
tained  In  a  causal  way.i.e.,  they  are  smoothing  residuals.  According  to  the 
ixed  representation,  some  of  the  r  pseudo  innovations  are  actually  passed 
rough  anti-causal  filters  instead  of  causal  ones  to  form  the  pseudo  processes, 
is  last  result  appears  closely  related  to  the  notions  of  direct  and  inverse  chan- 
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Hore  recently  developed  doubling  type  algorithms  [M3],  that  promise  to  be 

even  more  efficient,  especially  for  large  problems. 

Ve  present  several  results  that  allow  us  to  see  these  seemingly  different 
methods  as  part  of  a  whole  set  of  tools  that  are  closely  related.  The  new  results 
can  be  viewed  as  missing  links  to  a  larger  picture  that  is  gradually  emerging. 

First  we  present  a  new  Toeplitz  block  matrix  embedding  principle,  that  pro¬ 
vides  an  elegant  derivation  of  the  Generalized  Levinson  algorithm,  as  well  as  con¬ 
nections  to  other  problems,  such  as  the  inversion  of  Toeplitz  plus  Hankel 
matrices.  The  basic  idea  is  rather  simple;  it  is  well  known  that  the  inverse  of  a 
block  matrix  can  be  expressed  via  the  block  inversion  formulas  as  a  rational 
combination  of  the  component  blocks,  the  idea  is  to  reverse  that  process  and 
view  a  given  matrix  as  a  rational  combination  of  component  blocks,  where  the 
blocks  are  components  of  a  larger  matrix.  In  this  embedding  process  the  hope  is 
that  the  larger  matrix  displays  a  structure  that  is  simpler  than  the  original 
(smaller)  matrix. 

A  historical  example  of  this  type  is  one  of  Levinson's  attempts  to  solve  Toe- 
plitz  systems  in  [W].  By  partitioning  a  Toeplitz  matrix  into  a  two  by  two  Toeplitz 
block  matrix  and  applying  a  block  reduction  step,  he  winds  up  with  a  Toeplitz 
plus  Hankel  matrix  inversion  problem  of  half  size.  The  problem  was  abandoned 
at  this  stage.  Our  point  is  that  the  inverse  process  can  be  used  to  embed  the 
Toeplitz  plus  Hankel  matrix  inversion  problem  into  a  Toeplitz  block  matrix  or, 
via  an  interleaving  permutation,  into  a  block  Toeplitz  matrix.  Block  Toeplitz 
matrices  can  be  inverted  with  the  multivariate  Levinson-Whittle-Wiggins- 
Robinson  (LWR)  algorithm,  while  Toeplitz  block  matrices  are  inverted  using  the 
Toeplitz  block  version  of  the  LWR  algorithm,  reminiscent  of  two-dimensional  (2- 
D)  generalizations  of  the  LWR  algorithm,  see  e.g.  [LKM].  The  Toeplitz  block 
matrix  form  of  the  LWR  algorithm  can  be  extended  to  matrices  with  Toeplitz 
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as  auxiliary  quantities.  In  estimation,  these  quantities  are  interpreted  as  covari¬ 
ances  of  residuals  of  one  step  predictors.  In  the  doubling  algorithm  for  a- 
Toeplitz  matrices  [M3],  Schur  complements  of  higher  orders  have  to  be  con¬ 
sidered.  The  tree  associated  with  the  algorithm  is  binary  with,  sitting  at  the 

nodes,  matrices  of  order  1,  2,  4 . n  obtained  by  combined  partitioning  and 

Schur  complementation.  The  order  1  quantities  are  the  same  as  the  auxiliary 
quantities  (residuals)  associated  with  the  leftmost  (or  rightmost)  branch  of  the 
linear  tree.  The  order  2*  quantities  have  the  same  interpretation  taking  order 
2k  blocks  as  matrix  elements.  The  algorithm  is  function  recursive,  a  property 
that  appears  well  adapted  to  the  tree  machines  [BRO],  [SEQ].  The  a-Toeplitz 
matrix  is  inverted  using  a  depth-first  traversal  of  the  tree;  this  total  ordering 
implies  a  sequential  overall  computational  scheme,  although  every  (order) 
update  may  be  computed  in  a  parallel  fashion  (via  FFTs).  As  for  the  linear  tree, 
computations  are  further  reduced  if  some  quantities  at  the  same  level  in  the 
tree  are  known  a  priori  to  be  identical  or  have  low  rank  differences.  The 
significance  erf  such  a  structure  is  explored.  A  very  simple  example  is  given  by 
the  covariance  matrix  of  a  Brownian  motion;  then  only  one  branch  of  the  binary 
tree  has  to  be  considered.  For  possible  VLSI  architectures,  see  [ADM]. 

Mixed  and  Minimal  Representations  for  Toeplitz  and  Related  Systems  [DM1] 

Many  signal  processing,  estimation  and  control  problems  require  the  solu¬ 
tion  of  Toeplitz,  Hankel  or  related  systems  of  equations.  Such  solutions  can  be 
obtained  efficiently  via  the  use  of  several  existing  methods; 

The  Generalized  Levinson  algorithm,  see  e.g.  [FMKL],  for  matrices  close  to 
Toeplitz,  such  as  sums  of  products  of  Toeplitz  matrices. 

Ladder  or  Lattice  Form  realizations,  see  e.g.  [MLNV],  that  result  in  stable 
and  efficient  implementations. 


. 


trices  representative  ot  each  class.  If  the  matrix  is  Toeplitz,  all  contiguous  prin¬ 
cipal  submatrices  of  the  same  order  are  identical  so  that  only  one  branch  of  the 
tree  has  to  be  considered  and  the  computations  are  equivalent  to  the  Levinson 
algorithm  [LEV].  Since  order  updates  of  (distinct)  submatrices  of  the  same 
dimension  can  be  performed  simultaneously,  the  algorithm  is  well-suited  for 
parallel  processing.  More  precisely,  with  a  parallel  processing  implementation 
of  the  algorithm,  any  matrix  is  inverted  with  the  same  number  of  steps  as  a  Toe¬ 
plitz  matrix  of  the  same  size. 

In  many  important  examples  related  to  estimation  and  control,  contiguous 
principal  submatrices  of  the  same  order  have  a  low  rank  difference.  When  the 
difference  of  the  two  contiguous  principal  submatrices  of  order  n  —  1  has  (low) 
rank  a.  the  rank  of  the  difference  of  two  contiguous  principal  submatrices  of 
any  order  is  bounded  by  a  and  the  matrix  has  a  shift  low  rank  a  [Ml]  and  was 
referred  to  as  a-Toeplitz  in  [FMKL].  The  contiguous  principal  submatrices  of  an 
a-Toeplitz  matrix  are  all  distinct  in  general,  hence,  since  all  the  branches  of  the 
tree  would  have  to  be  explored,  the  algorithm  developed  in  [DGKl]  would  not 
invert  efficiently  such  a  matrix.  However,  since  two  contiguous  principal  subma¬ 
trices  differ  by  a  low  rank  correction,  (a  representation  of)  the  inverse  of  one  is 
easily  deduced  from  the  inverse  of  the  other;  this  is  a  one  step  time  update. 
Thus  using  both  "vertical"  (order)  and  "horizontal"  (time)  updates,  only  two 
principal  submatrices  of  each  order,  i.e.,  two  branches  of  the  tree,  need  to  be 
considered.  This  provides  a  nice  interpretation  of  the  Generalized  Levinson  algo¬ 
rithm  [FMKL]  [FKML].  It  is  worth  mentioning  at  this  point  that  the  Generalized 
Levinson  algorithm  for  inverting  a-Toeplitz  matrices  may  also  be  derived  from 
the  (block)  Levinson  algorithm  applied  to  a  particular  (a  +  l)  x  (a  +  1)  block 
Toeplitz  matrix  [DMl]. 

The  algorithms  considered  up  to  now  involve  first  order  Schur  complements 
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4.3  Nearly  Stationary  or  finite  Rank  Processes 

In  this  section  we  present  summaries  of  our  continuing  work  on  nearly  sta¬ 
tionary  processes.  We  made  significant  advances  towards  a  better  understand¬ 
ing  of  such  processes  and  the  algorithms  related  to  their  estimation.  Embed¬ 
ding  techniques  play  a  central  role  in  our  treatment  of  nearly  stationary  pro¬ 
cess.  We  were  able  to  show  that  a  large  class  of  problems  involving  non¬ 
stationary  processes  of  a  certain  type,  can  be  embedded  in  a  multichannel  sta¬ 
tionary  problem  The  papers  summarized  below  provide  a  comprehensive 
framework  for  deriving  efficient  estimation  algorithms  for  these  processes. 

Significant  progress  was  also  made  towards  achieving  a  stochastic  interpretation 
of  the  generalized  Levinson  algorithm.  We  have  a  much  better  understanding  of 
this  algorithm  and  of  its  different  implementations.  It  is  clear  now  that  the  time 
and  order  update  equations  of  the  generalized  Levinson  algorithm  lend  them¬ 
selves  naturally  to  distributed  implementation.  A  more  detailed  discussion  of 
this  point  is  presented  next. 

A  Tree  Classification  of  Algorithms  for  Toeplitz  and  Related  Equations  [MD1] 

The  development  of  VLSI  realizations  of  tree  machines,  see  e.g.,  [BRO], 

[SEQ]  motivates  a  study  of  the  structure  of  fast  matrix  inversion  algorithms 
from  a  tree  viewpoint,  see  also  [ADli]. 

An  algorithm  for  efficiently  inverting  matrices  "close  to  Toeplitz"  has  been 
recently  proposed  in  [DGK1].  In  that  context,  an  n  x  n  matrix  is  considered 
close  to  Toeplitz  if  its  contiguous  principal  submatrices  of  any  given  order 
belong  to  only  a  few  classes,  where  two  matrices  belong  to  the  same  class  if  they 
are  either  identical  or  the  reverse  of  each  other.  The  tree  underlying  the  algo¬ 
rithm  is  the  linear  tree  associated  with  the  contiguous  principal  submatrices  of 
order  1  to  a  The  algorithm  performs  (one  step)  order  updates  to  invert  subma- 
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4.2  Two-Dimensional  Systems 

Our  earlier  studies  in  2-D  systems  have  revealed  that  it  is  very  difficult  to 
generalize  most  of  the  essential  results  for  1-D  systems,  such  as  signal  process¬ 
ing  techniques  associated  with  only  time  dependent  signals.  This  led  us  to  reex¬ 
amine  the  1-D  techniques  and,  in  particular,  study  very  carefully  the  estimation 
problems  for  processes  with  finite  displacement  rank,  see  e.g.,  [DMl]-[DGMV], 
since  we  realized  that  this  theory  can  be  extended  to  tackle  the  statistical  filter¬ 
ing  of  non-stationary  random  fields. 

The  work  of  Livshits  on  operator  colligations  [L]  and,  more  specifically,  on 
vector  colligations  [LY]  as  well  as  some  more  recent  papers  [L2],  [KRA]  lay  the 
ground  to  such  a  generalization.  The  theory  is  based  on  the  use  of  data  col¬ 
lected  along  paths  and  is  in  that  sense  akin  to  an  approach  of  Willsky  et  al.  [WS]. 
As  a  consequence  it  can  be  used  for  estimation  not  only  for  random  fields  on  the 
plane  or  any  Euclidean  space  but  also  for  instance  on  the  sphere;  hence  is 
naturally  adapted  to  geo-location  problems. 

Although  Livshits  does  not  study  the  estimation  problem,  he  introduces  the 
basic  concepts  necessary  for  filtering  based  on  covariance  information  only. 
Filtering  based  on  second-order  statistics  is  very  appealing  from  a  practical 
standpoint  since  usually  estimates  of  higher  order  moments  are  not  reliable. 
Furthermore  local  state-space  models,  such  as  Roesser’s  for  2-D  random  fields, 
see  e.g.,  [KLUK],  do  not  need  to  be  introduced.  The  notions  of  reflection 
coefficients  and  ladder  structures  for  finite  rank  random  fields  appear  as  a  key 
concept  from  both  theoretical  and  computational  viewpoints.  An  important 
issue  is  the  determination  of  the  relationship  between  the  reflection  coefficients 
obtained  from  this  approach  and  the  ones  introduced  by  Marzetta  [MAR]  and 
Delsarte  et  al.  [DGK2]-[DGK3]  in  the  particular  case  of  stationary  2-D  random 
fields  with  data  collection  using  line-by-line  scanning. 
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dimensional  systems. 

We  are  currently  working  on  developing  basic  results  for  multi-dimensional 
systems  which  will  provide  a  general  framework  for  handling  DSN-type  problems. 
Preliminary  research  has  indicated  that  there  are  non-trivial  differences 
between  the  2-D  and  3-D  cases,  which  makes  direct  extensions  of  known  results 
more  difficult. 

Another  fundamental  issue  that  needs  to  be  addressed  is  the  assumption  of 
stationarity,  which  is  a  basis  for  practically  all  current  estimation  techniques. 
The  time  varying  nature  of  acoustic  propagation  models  and  the  various  events 
that  can  occur  in  a  DSN  system  are  highly  non-stationary.  It  is  imperative, 
therefore,  to  develop  estimation  techniques  capable  of  handling  non-stationary 
processes.  Using  the  notion  of  nearly  stationarity,  we  obtained  a  useful  charac¬ 
terization  of  non-stationary  processes,  that  lends  itself  to  analysis  and  algorithm 
development  Significant  advances  have  been  made  in  our  understanding  of 
non-stationarity  and  a  number  of  interesting  results  were  obtained.  Another 
important  aspect  of  our  work  on  non-stationary  processes  iB  related  to  the  pos¬ 
sibility  of  their  implementation  by  distributed  algorithms.  Our  analytical  frame¬ 
work  turns  out  to  provide  a  natural  problem  decomposition,  that  can  be  utilized 
for  distributing  the  required  computations  among  multiple  processors. 


(phase  shifts,  FFT-s)  fall  into  this  category. 

The  special  structure  of  the  source  location  problem,  in  a  homogeneous 
propagation  media,  and  the  assumption  of  Gaussian  statistics  makes  possible 
the  separation  of  spatial  and  temporal  estimation.  However,  this  separation  fails 
for  more  complicated  propagation  models.  Even  in  relatively  simple  cases,  this 
type  of  approach  is  not  completely  satisfactory  where  multiple  sources  are 
present. 

Our  approach  has  therefore  been  to  consider  a  multi-dimensional  estima¬ 
tion  problem.  In  fact,  when  only  a  finite  number  of  discrete  sensors  are  used, 
the  problem  is  multi-channel,  rather  than  multi-dimensional. 

The  system  consisting  of  several  sources  and  sensors  can  be  viewed  as  a 
multi-input  multi-output  (MIMO)  system.  Under  certain  assumptions  on  the 
source  spectra  and  the  propagation  model,  the  vector  time-series  observed  by 
the  sensors  can  be  modeled  as  a  multi-channel  ARMA  process.  By  using  various 
least  squares  or  maximum  likelihood  parameter  estimation  techniques,  it  is  pos¬ 
sible  to  fit  a  set  of  ARMA  coefficients  to  the  measured  data  We  have  been  able  to 
show  that  source  locations  and  their  spectra  can  be  obtained  relatively  easily 
from  these  coefficients.  In  general,  there  exists  a  mapping  between  the  proper¬ 
ties  of  the  source  location  problem  into  the  structure  of  the  ARMA  coefficients. 
The  details  of  this  mapping  and  the  structure  it  imposes  on  the  MIMO  system  is 
presently  under  investigation. 

As  mentioned  above,  multichannel  estimation  is  a  special  case  of  multi¬ 
dimensional  estimation.  A  large  body  of  results  is  currently  available  for  2-D  sig¬ 
nal  processing.  These  results  were  developed  mainly  in  the  context  of  image 
processing,  where  two  spatial  dimensions  are  considered.  Relatively  little  seems 
to  have  been  done  in  the  case  where  one  spatial  dimension  and  one  temporal 
dimension  are  involved.  Even  less  seems  to  be  known  about  estimation  in  higher 
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4.  BASIC  RESEARCH  ON  DISTRIBUTED  ALGORITHMS 

4.1  Orel  lien 

The  analysis  and  design  of  a  distributed  sensor  network  involves  a  variety  of 
issues.  Very  little  basic  research  has  been  done  on  the  fundamental  mathemati¬ 
cal  problems  that  need  to  be  solved  in  order  to  provide  a  systematic  approach 
to  DSN  system  design.  HVhile  most  of  our  work  concentrated  on  the  more 
immediate  problem  of  developing  signal  processing  tools  for  the  DSN,  we 
devoted  part  of  our  effort  to  more  basic  research.  In  this  section  we  discuss  two 
research  areas  of  this  kind:  multi-dimensional  estimation,  and  nearly  stationary 
processes. 

Multiple  sensor  measurements  can  be  considered  as  samples  of  a  time- 
space  function  y{t,£\p)  generated  by  the  source,  where  t  represents  the 
sensor  location  in  3-D  space,  and  p  is  a  vector  of  source  parameters.  This  func¬ 
tion  may  represent  electromagnetic  waves  reflected  from  an  object  or  acoustic 
waves  generated  by  aircraft  engines.  The  parameter  vector  p  may  include 
source  location,  velocity,  and  spectral  characteristics.  Thus  the  problem  is  to 
estimate  the  unknown  parameter  vector  p  from  a  set  of  noisy  measurements. 

Current  approaches  to  this  problem  are  almost  exclusively  based  on  reduc¬ 
ing  the  multi-dimensional  estimation  problem  to  a  sequence  of  one  dimensional 
problems.  Typical  array  processing  techniques  involve  forming  multiple  beams 
by  performing  linear  operations  on  the  sensor  outputs.  This  linear  operation 
reduces  the  dimension  of  the  problem  by  effectively  concentrating  on  a  single 
spatial  direction  at  a  time.  The  ensuing  signal  processing  involves  a  scalar  time 
series  and  thus  standard  estimation  and  filtering  techniques  can  be  applied.  All 
the  different  versions  of  beam-forming  and  TDOA  estimation  techniques,  whether 
implemented  in  the  time  domain  (delays,  correlators)  or  in  the  spectral  domain 
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Efficient  Solution  of  Lyapunov  Equations  [Pll  J.  [HH] 

The  Levinson  algorithm  for  fitting  an  autoregressive  (AR)  model  to  a  given 
covariance  matrix,  has  many  nice  properties.  In  addition  to  being  computation¬ 
ally  efficient  is  provides  not  only  a  set  of  AR  coefficients  but  also  the  related 
reflection  coefficients  and  the  backwards  predictor  coefficients.  In  fact,  the 
Levinson  algorithm  yield  two  Cholesky  factorizations  of  the  inverse  covariance 
matrix.  In  many  applications  we  encounter  the  inverse  problem:  given  a  set  of 
AR  coefficients,  find  the  covariance  sequence,  or  the  reflection  coefficients.  In 
DSN  applications  the  problem  arises  in  the  context  of  translating  "standard" 
parameters  (i.e.,  parameters  of  models  given  in  direct  form)  into  the  ladder 
parameters  required  by  our  algorithms.  A  similar  problem  is  encountered  when 
transforming  left  to  right  MfD's. 

The  solution  of  the  inverse  problem  is  relatively  easy  in  the  scalar  case. 
The  reason  is  that  the  backward  predictor  is  then  just  the  reversed  order  for¬ 
ward  predictor.  Therefore,  it  is  possible  to  run  the  Levinson  algorithm  back¬ 
wards  starting  from  the  given  full-order  solution  and  going  back  to  zero  order. 
Unfortunately,  this  procedure  does  not  carry  through  to  the  matrix  case,  since 
the  backward  predictor  is  not  longer  a  simple  function  of  the  forward  predictor 
coefficients.  Finding  the  covariance  sequence  requires  solution  of  the  discrete 
Lyapunov  equation.  Straightforward  solutions  of  this  equation  are  computation¬ 
ally  expensive  for  high  order  systems.  In  [PM]  and  [HM]  we  develop  novel  solu¬ 
tion  techniques  that  utilize  the  problem  structure  to  significantly  reduce  com¬ 
putational  requirements.  Three  efficient  algorithms  sure  presented.  One  of  the 
methods  is  based  on  a  recently  described  method  of  inverting  matrices  that  are 
sums  of  block-Toeplitz  and  block-Hankel  matrices  [FM3].  The  procedures  are 
then  shown  to  yield  a  stability  test  for  the  given  autoregressive  model. 
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By  properly  choosing  the  delay  D  (or  the  pre-filter),  the  wide  band  components 
at  the  two  inputs  become  decorrelated  while  the  narrow  band  signal  remain 
correlated.  Since  the  coefficients  of  the  filter  are  updated  to  minimize  the  error 
power,  the  predictor  and  error  outputs  can  be  applied  to  line  enhancement  and 
signal  whitening  respectively.  Most  analyses  of  the  ALE  were  for  a  white  noise 
background,  while  very  few  results  exist  for  the  colored  noise  and  the  whitening 
performance.  Satorius  and  Zeidler  considered  an  autoregressive  moving  aver¬ 
age  (ARMA)  input  process  [SZ]. 

In  [NM],  a  matrix  formulation  was  used  to  derive  the  optimal  coefficients 
and  frequency  response  of  the  ALE  and  exact  least-squares  predictors,  for  inputs 
of  real  or  complex  sinusoids  in  general  additive  colored  noise.  The  derivation  of 
the  amplitude  gains  of  the  output  sinusoids  generalized  known  results  for  the 
white  noise  case.  In  low-pass  and  high-pass  background  noise  the  amplitude 
gains  become  essentially  monotonically  increasing  and  decreasing  functions 
respectively  of  the  sinusoid  frequencies.  The  amplitude  distortion  that  is  intro¬ 
duced  by  the  correlated  noise  can  be  reduced  by  choosing  a  larger  number  of 
coefficients,  L.  In  evaluation  of  the  performance  of  the  predictors  in  whitening 
applications,  an  upper  bound  on  the  output  SNR  was  found  when  filtering  a  white 
signal  that  had  been  correlated  by  multiple  stationary  sinusoids. 

To  enable  filtering  of  nonstationary  complex  and  multichannel  data,  a  com¬ 
plex  vector  version  of  the  ladder  algorithm  is  developed.  It  can  be  used  to 
implement  the  complex  ALE,  as  well  as  the  noise  canceller  or  noise  inversion.  A 
major  advantages  of  this  algorithm  is  its  rapid  rate  of  convergence  to  the  exact 
least-squares  solution,  hence  it  can  actually  achieve  the  upper  bounds  of  the 
SNR  performance. 


used:  some  fraction  of  the  alpha  pseudo  innovations  is  actually  passed  through 
anti-causal  or  adjoint  filters  instead  of  causal  ones.  This  mixed  representation 
avoids  a  particular  signature  problem.  Physically,  the  necessary  algebra  is  very 
much  reminiscent  of  modern  physics;  the  different  pseudo  processes  have  simi¬ 
lar  properties  as  particles  and  anti-particles,  in  the  sense  that  they  can  be  con¬ 
verted  from  real  causal /anti-causal  representations  to  imaginary  anti- 
causal/caused  representations.  Other  evidence,  such  as  the  fact  that  Livshits 
definition  of  the  (displacement)  rank  of  a  process  generalizes  to  multi¬ 
dimensions,  makes  the  physical  insight  even  more  interesting.  The  implications 
of  these  results  to  continuous  time  models  and  the  solution  of  Fredholm  integral 
equations  are  outlined. 

As  a  last  topic  we  mention  the  direct  identification  of  alpha  stationary 
processes.  Because  of  the  multiplicity  issue  the  question  of  degree  of  freedom  of 
models  for  such  processes  has  to  be  raised.  A  direct  answer  would  probably  be 
quite  difficult;  however,  a  different  parametrization  of  the  models  using  so-called 
ladder  or  lattice  form  realizations,  see  e.g.  [MLNV],  [ML2],  [LEE],  lead  to  a  very 
promising  alternative.  We  only  note  here  that  the  gaussian  log-likelihood  func¬ 
tions  and  their  derivatives  with  respect  to  the  model  parameters  are  simple 
functions  that  can  be  straight  forwardly  analysed. 

Distortion  Measures  of  Finite  Rank  Processes  Via  Ladder  Forms  [LM4] 

It  has  been  known  for  some  time  that  distortion  measures  of  second  order 
processes  can  be  evaluated  using  state-estimation  and  related  techniques. 
Measures  such  as  the  Bhattacharyya,  the  Kullback-Leibler,  and  the  related 
Itakura-Saito  for  second  order  processes  involve  determinants,  traces  and  ratios 
of  covariance  matrices,  see  e.g..  [KA3].  In  many  applications,  it  is  of  interest  to 
compute  these  measures  efficiently.  In  speech  processing,  these  measures  are 
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important,  for  instance,  in  the  design  of  low  rate  speech  encoder  [MBG],  as  well 
as  for  speech  recognition  both  for  the  discrete  utterance  [1TA],  and  continuous 
speech  cases  [JEL],  In  pattern  recognition,  these  measures  are  important  in 
feature  extraction  and  clustering  analysis.  In  digital  communication,  these 
measures  arise  in  channel  capacity  calculations. 

In  many  applications,  one  may  be  interested  in  efficient  parametrizations  of 
such  measures,  either  to  simplify  the  calculations  and  design,  or  in  the 
hardware  implementation.  Efficiency  is  of  interest  if  the  optimization  of  such 
measures  has  to  be  done  "on-line"  or  "adaptively".  Such  an  optimization  has 
become  feasible,  due  to  the  availability  of  very  large  scale  integrated  circuits 
(VLSI). 

It  is  well  known  that  distortion  measures  of  stationary  processes  involve 
covariance  matrices  and  kernels  of  the  Toeplitz  or  displacement  type.  These 
measures  require  either  the  solution  of  linear  equations  involving  such  opera¬ 
tors,  or  the  computation  of  triangular  factors  or  volterra  kernels.  They  can  be 
obtained  efficiently  via  the  use  of  existing  methods,  such  as  the  Krein-Le vinson 
type  algorithms. 

The  first  aim  of  this  paper  is  to  emphasize  that  the  stationarity  is  not 
necessary  tor  obtaining  computational  benefits  and  ease  of  parametrization; 
that  is,  we  are  now  able  to  efficiently  compute  distortion  measures  for  processes 
that  deviate  from  stationarity  in  a  certain  sense.  An  index  of  stationarity,  called 
"shifter  low  rank"  (displaced)  rank  of  alpha,  has  recently  been  introduced, 
[Mo].[FMKL]  [KLM]  [KKM],  with  the  significance  that  the  larger  the  value  of  alpha 
the  more  non-stationary  a  process  is.  In  addition,  the  complexity  of  process 
models  and  calculations  involving  such  processes  increases  only  proportional  to 
this  index.  In  the  discrete  case,  alpha  can  be  used  to  measure  the  distance  of  a 
matrix,  e.g.,  a  covariance  matrix,  from  being  Toeplitz.  Our  earlier  results  showed 
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that  Kre in-Levinson  and  other  algorithms  for  Toeplitz  operators  can  be  general¬ 
ized  to  such  alpha-stationary  processes,  as  well  as  a  new,  and  even  more 
efficient  doubling  type  algorithms  [M3]. 

Secondly,  we  show  that  ladder  canonical  forms,  which  are  realizations  of 
processes  based  on  partial-correlation  coefficients,  have  many  advantages  in 
computing  distortion  measures.  These  ladder  forms  in  fact  represent  a 
"natural"  canonical  parametrization  of  many  of  the  measures,  in  the  sense  that 
the  measures  become  simple  functions  of  the  ladder  form  parameters.  They 
typically  involve  weighted  averages,  including  logarithms.  For  instance  the 
gaussian  log-likelihood  derivative  with  respect  to  these  ladder  parameters  is  a 
linear  function  in  these  parameters!  Moreover,  the  ladder  forms  have  many 
other  nice  features,  such  as  stability  and  automatic  seeding  properties. 

Several  examples  involving  distortion  measures  are  given,  and  their  physi¬ 
cal  significance  and  applications  to  communication  and  estimation,  e.g.  speech 
processing  and  digital  communication  is  discussed. 

Efficient  Inversion  Formulas  far  Sums  of  Products  of  Toeplitz  and  Hankel 
Matrices  [FH3] 

The  problem  of  solving  linear  equations,  or  equivalently  of  inverting 
matrices,  arises  in  many  fields.  Many  techniques  for  doing  this  have  been 
developed,  but  all  of  them  have  the  characteristic  that  it  takes  in  the  order  of 
N 9  operations  (multiplications  and  additions)  to  invert  a  general  N  x  N  matrix. 
A  great  deal  of  work  has  been  devoted  to  finding  more  efficient  ways  of  inverting 
matrices,  by  utilizing  any  special  structure  that  they  may  posses.  This  special 
structure  can  take  various  forms,  but  there  are  many  problems  in  mathematical 
physics  and  statistics  in  which  the  matrices  can  be  shown  to  be  Toeplitz  or  block 
Toeplitz,  i.e.. 
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R  =  ['li-i)]  ■  0  *  i .}  <  N  ,  (4.1) 

where  r  is  a  p  x  p  block  entry.  Some  examples  for  applications  involving  Toe- 
plitz  matrices  are  least  squares  estimation  problems  involving  convolutions  and 
integral  equations  with  a  difference  kernel.  Hankel  matrices  represent  another 
type  of  special  structure,  characterized  by 

R  =  [r(t+i)]  •  Os z  ij  *  N  .  (4.2) 


Several  recursive  algorithms  have  been  developed  for  the  inversion  of  Toe- 
plitz  and  Hankel  matrices,  requiring  only  0(NV)  arithmetic  operations.  The 
extension  of  this  work  to  block-Toeplitz  and  block-Hankel  matrices  of  size  Np  x 
Np  (Le.,  NxN  blocks  of  size  p  x  p),  shows  that  the  inverse  can  be  obtained  in 
Oip^N*).  Recently  some  new  algorithms  have  been  developed  that  compute  the 
inverse  of  a  Toeplitz  matrix  even  faster,  in  0(NlogzN)  operations. 

While  Toeplitz  and  Hankel  matrices  appear  in  numerous  applications,  it 
often  happens  that  practical  problems  involve  matrices  with  a  more  complex 
structure.  As  an  example,  the  inverse  of  a  Toeplitz  matrix  is  no  longer  Toeplitz 
but  it  retains  a  very  simple  structure.  In  previous  work  we  have  shown  that 
efficient  inversion  formulas  can  be  obtained  for  matrices  that  are  "close  to  Toe¬ 
plitz''  in  some  sense.  More  precisely,  we  considered  the  class  of  matrices 
represented  by 


where 


T°  .=  T  +  t  OiUUi  , 

t»i 


T  =  full  Toeplitz  matrix 


(4.3) 


It(C4)  =  lower  (upper)  triangular  Toeplitz  matrix, 


ff4  =  ±  1 
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and  a  is  the  smallest  integer  for  which  this  decomposition  holds.  We  have 
shown  that  the  number  a  associated  with  a  given  matrix  is  unique,  and  given  by 

a  =  rank  d[r*]  (4.4) 


R  -  [rij]  .  0  *  i.j  <  N  . 


(4.5b) 


We  will  refer  to  a  matrix  of  this  type  as  having  a  shift-low  rank  or  displace¬ 
ment  rank  of  a,  or  more  briefly,  as  an  a-Toeplitz  matrix.  This  class  of  matrices 
has  a  nice  closure  property  in  the  sense  that  the  inverse  of  a-Toeplitz  matrix  is 
also  a-Toeplitz. 

The  natural  extension  of  these  results  to  matrices  that  are  "close  to 
Hankel,"  involves  the  following  type  of  representation: 

Ha  =  Tl  +  £  Li  Ui  7  =  T*I  =  T  ■  .  (4.6) 

i*i 

where  T.  Z.  U,  7t  are  defined  in  the  same  way  as  T,  L,  U,  o.  In  other  words, 
the  a-Hankel  matrix  is  simply  a  column-permuted  version  of  an  a-Toeplitz 
matrix.  Thus,  a  trivial  way  of  inverting  a  a-Hankel  matrix  is  by  transforming  it 
into  an  a-Toeplitz  matrix  and  then  applying  the  algorithm  given  for  instance  in 
[FMKL],  i.e., 

+  t  *  A  d’  7  •  (4.7) 

i*l 


However,  it  is  also  possible  to  invert  the  a-Hankel  matrix  directly.  The 
difference  between  these  two  types  of  inversion  formulas  lies  in  the  sequence  of 
submatrices  that  are  recursively  inverted.  This  subtle  difference  makes  the 
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trivial  solution  (4.7)  inapplicable  in  many  situations,  e.g..  if  H*  =  /T*  ©  Hl. 

In  [FM3]  we  extend  the  class  of  matrices  for  which  efficient  inversion  is  pos¬ 
sible,  to  include  mixtures  of  Toeplitz  and  Hankel  matrices.  We  consider 
matrices  of  the  form 

R  -  T  4  2  Oi  Li  Vi  4  T  T  4  £  Li  Vt  T  =  T*  +  H*  ,  (4.8) 

(■i  t«i 

Le.,  sums  of  a-Toeplitz  and  0-Toeplitz  matrices.  This  class  of  matrices  has  a  clo¬ 
sure  property  in  the  sense  that  R~l  has  a  representation  of  the  type  given  in 
Equation  (4.8). 

The  basic  idea  underlying  our  results,  is  to  embed  the  a-Toeplitz //?-Hankel 
matrix  in  a  larger  matrix  that  has  a  structure  for  which  efficient  inversion  for¬ 
mulas  have  already  been  developed.  Two  such  embeddings  are  explored:  We 
first  show  that  an  a-Toe plitz / /9-Hanke l  matrix  R  with  p  x  p  entries  (i.e., 
R  -  [r4j],  r<j  is  a  p  x  p  matrix),  can  be  embedded  in  an  ( a +/?) -Toeplitz 
matrix  with  2p  x  2p  blocks.  Applying  the  inversion  formulas  derived  in  [FMKL] 
to  this  matrix  results  in  an  algorithm  which  requires  0(4(a  +  jJ  +2)pz7Vz) 
operations.  Second  we  present  a  way  of  embedding  the  a-Toepli tz //3-Hanke  1 
matrix  in  a  block  Toeplitz  matrix  with  2(a  4  p  +2)px2(a  4  p  2 )p  entries. 
Applying  the  multichannel  Levinson  algorithm  to  this  matrix  requires 
0(B(a  4  p  42)3p3JV2)  operations.  However,  by  using  the  special  structure  of  the 
block  entries  the  algorithm  can  be  further  simplified,  resulting  in  an  algorithm 
equivalent  to  the  one  mentioned  before.  Further  reductions  in  computational 
requirements  can  be  obtained  by  applying  the  class  of  doubling  formulas  (i.e., 
the  0(moi*N)  algorithms  described  in  [GY]  [AHV]  [M3]  to  the  embedded  block 
Toeplitz  matrix. 
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&  DEVELOPMENT  of  SIGNAL  PROCESSING  FACILITIES 

9|nat  Processing  Workbench 

An  essential  part  in  the  process  of  developing  new  techniques  for  source 
parameter  estimation  is  the  testing  of  proposed  algorithms  by  computer  simula¬ 
tion.  Simulations  allow  us  to  study  various  issues  related  to  the  performance  of 
the  algorithms  and  to  gain  insights  into  their  behavior.  Analysis  alone  is  not 
adequate,  especially  when  recursive  stochastic  algorithms  are  concerned. 
Simulations  are  an  invaluable  tool  for  studying  issues  such  as  convergence  rates, 
robustness,  and  estimation  accuracy.  Therefore,  establishment  of  a  signal  pro¬ 
cessing  workbench  that  will  allow  easy  comparison  and  modification  of  signal 
processing  steps  and  input  data  sequences  in  important.  In  particular,  this  will 
facilitate  the  exchange  of  processing  tools  and  data  bases  among  DSN  research¬ 
ers,  and  maximize  the  impact  of  our  research  results  on  the  design  of  a  DSN  sys¬ 
tem. 

In  the  proposed  workbench,  the  processing  modules  function  in  a  stand¬ 
alone  fashion  with  standard  input/output  formats  such  that  complex  processing 
functions  are  built  by  the  simple  cascading  of  modules.  In  the  UNIX  operating 
system,  precompiled  stand  alone  modules  can  be  connected  together  using  a 
system  command  called  a  pipe.  Thus  signal  processing  can  be  changed  quickly 
without  recompiling  or  linking,  just  forking  to  the  appropriate  module.  The  sig¬ 
nal  processing  workbench  should  readily  allow  inputs  derived  from  actual  exper¬ 
iments  or  simulated  inputs  and  allow  for  their  perturbation  by  various  types  of 
noise. 

Observing  and  modifying  the  internal  parameters  of  a  signal  processing  pro¬ 
cedure  while  data  is  being  processed  greatly  enhances  the  debugging  and  under¬ 
standing  of  signal  processing  algorithms.  The  version  7  UNIX  operating  system 
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for  the  VAX  computer  incorporates  a  debugger  that  allows  the  internal  parame¬ 
ters  of  a  signal  processing  submodule  to  be  observed  and  modified  as  the  data  is 
being  processed.  The  current  UNIX  debugger  applies  only  to  procedures  within  a 
main  program.  For  cascaded  stand  alone  programs,  this  facility  must  be 
extended.  Some  experience  with  the  problem  of  adding  this  capability  has  been 
obtained  on  the  DEC-11/34  through  the  development  of  our  own  debugger.  This 
experience  will  be  applied  to  extensions  of  the  debugger  for  the  VAX. 

In  several  theoretical  and  applied  problems  we  are  running  up  against  a 
boundary  characterized  by  the  need  of  manipulating  formulas.  This  raises  the 
urgent  need  for  obtaining  software  to  handle  algebraic  formulas  symbolically. 
We  made  an  effort  to  obtain  a  copy  of  MACSYMA  for  this  purpose;  however,  we 
were  unsuccessful  so  far,  due  to  lack  of  cooperation  from  the  originators  of 
MACSY14A. 

DSN  Computer  Facility 

The  recently  expanded  computer  facility  for  investigating  DSN  strategies  is 
centered  around  a  Digital  Equipment  Corporation  VAX  11/780,  see  Figure  5.1, 
running  under  the  UNIX  operating  system.  The  installed  operating  system  is  the 
latest  available  from  the  Electrical  Engineering  and  Computer  Science  Depart¬ 
ment  of  the  University  of  California  at  Berkeley,  the  Fourth  Berkeley  Software 
Distribution  (4BSC)  of  UNIX  version  7,  released  November  1960.  Our  current  VAX 
computer  hardware  incorporates  a  2.5  Megabyte  main  memory  subsystem  and 
80  Megabytes  of  disk  storage,  provided  by  a  single  RM03  disk  drive.  This  storage 
is  to  be  augmented  with  an  Ampex  300  MB  disk  drive  by  early  1961.  Our  com¬ 
puter  is  interconnected  with  a  DEC  11/70  (SU-ISL)  and  a  VAX  11/780  (SU-Mount 
St.  Helens)  and  to  the  ARPANET.  The  SU-DSN  VAX  became  operational  in 
December  1980  and  the  SU-ISL  11/70  was  brought  up  in  November  1980,  replac- 
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ng  a  much  overworked  DEC  11/34  system.  The  11/70  is  running  UNIX  version  6. 
installed  on  the  VAX  is  the  new  RAND  text  page  editor  (called  E)  which  gives  the 
jser  powerful  functions  for  moving  text  within  a  file  and  provides  greatly 
mproved  recovery  of  working  data  sets  after  system  crashes. 

In  the  last  two  years  prior  to  the  acquisition  of  the  VAX  and  11/70,  several 
somputers  were  used  to  test  algorithms  and  to  evaluate  their  performance 
under  simulated  conditions.  Among  these  facilities  used  were  the  ISL  11/34,  the 
VAX  (Mt.  St.  Helens)  which  became  operational  on  April  19B0,  and  the  facilities  of 
the  Stanford  Artificial  Intelligence  Laboratory. 

The  additional  300  MByte  disk  drive,  to  be  delivered  in  the  spring  of  1961,  is 
necessary  to  handle  the  quantity  of  data  expected  in  the  typical  DSN  scenario. 
Further  establishment  of  an  interactive  graphics  system  will  enhance  our  work 
on  application  of  image  reconstruction  techniques  to  the  DSN  problem.  Two 
high  resolution  bit  map  display  systems  are  being  added  to  our  VAX.  one  can  be 
upgraded  to  gray  scale  and  color  later.  We  have  ordered  a  laser  printer  to 
upgrade  the  quality  of  text  and  hardcopy  graphics  output.  In  order  to  efficiently 
perform  image  reconstruction  and  target  tracking  using  data  from  multiple  sen¬ 
sor  nodes,  our  fast  access  memory  should  be  extended. 


-87- 


a)  SU-DSN: 


VAX  11/780 
UNIX  V7 
4-th  BSD 


LS^llJ - 


memory 

subsystem 

disk 

2.5 

RM03 

Mb 

— « — 

80  Mb 

TU  77 

CPU 

FPA 

cache 

tape  drive 
800/1600  bp 
125  ips 


massbus 

adapters 


uni bus  C  UBA 
adapter  j 


Jenson/Variar 

jlotter 


ARPANET 
9  SUMEX 


DL-11 


Ampex 
disk  9300 
storage  300  Mb 


controller 


<■ 


terminals  &  modems  ^  _ 

(SU-ISL)  DEC  11/70  _<rl —  1  ■  : 

(SU-VLSI)  VAX  11/780  = - = 

SU-ISL  FACILITIES  DZ_n  (5) 

DEC  11/70  375kb  CPU  memory  &  FPU-A 
2  Ampex  9300  300Mb  disks 
Digldata  tape  drive  800/1600  bpi  45  Ips 
Printronix  line  printer  3001pm 
2  Diablo  1640  terminals 
Tektronix  4013  graphics  terminal 
dial  out  modems 

terminal  Interconnect  facility  to  SU-DSN,  SU-ISL,  su-MT.  ST.  HELENS  and  other 

*  micro-computer  systems 


Figure  5.1 


SU-DSN  COMPUTER  FACILITIES 


A*'  -.’.i'-AV.  j  .'l.l'*,.  I'J  •l.'TJ.'V  21 . 


i.i.i.  r.'  i 


-88- 

nn  and  Hardware  Programming  Efforts 

3ur  early  efforts  in  software  development  were  aimed  at  establishing  a  pro- 
i  library  of  system  identification  and  parameter  estimation  procedures. 
Library  includes  standard  recursive  identification  methods,  prewindowed 
sr  algorithms  in  both  AR  and  ARMA  forms,  and  newly  developed  covariance 
sliding  window  ladder  form  algorithms.  Normalized  and  unnormalized  ver- 
i  of  the  algorithms  have  been  written.  The  ladder  form  algorithms  were 
en  primarily  in  "C"  while  the  other  routines  were  written  in  Fortran  or 
[SAIL,  the  language  of  the  Stanford  Artificial  Intelligence  Laboratory. 

A  speech  modeling  system  was  previously  developed  around  the  prewin- 
sd  ladder  form  as  described  in  [MLl].  The  reflection  coefficients  which 
meterize  the  speech  spectrum  were  determined  recursively  with  updates 
y  new  speech  sample.  At  the  desired  transmission  rate,  the  best  single 
e  of  a  reflection  coefficient  is  chosen  The  likelihood  parameter  from  the 
ration  coefficient  updating  equation  was  used  to  identify  the  pitch  pulse  loca- 
.  Event  detection  algorithms  that  indicate  a  sudden  change  in  the  underly- 
iignal  structure  evolved  from  the  pitch  extraction  technique.  ‘ 

The  prewindowed  ladder  form  has  also  been  the  basis  of  an  investigation 
efficient  hardware  implementations  using  the  C0RD1C  computation  tecb- 
e  to  reduce  the  complexity  of  the  algorithm.  The  study  led  to  a  series  of 
I  chip  designs,  that  turned  out  to  be  extremely  promising  for  general  real 
i  digital  signal  processing.  A  data  modem  using  the  ladder  form  was  also 
nined. 

The  linear  equation  technique  for  solving  the  source  location  problem  from 
t  difference  of  arrival  data  has  been  verified  by  computer  simulation.  The 
:ts  of  imprecise  sensor  location  and  noisy  TDOA  information  have  also  been 
dated. 
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Work  on  modification  of  medical  image  reconstruction  techniques  for  the 
DSN  scenario  is  continuing.  The  minimum  variance  approach  to  Image  Recon¬ 
struction  developed  by  Wood  [WOO],  Fortes  [FOR]  and  Nunes  [NUN]  originally  was 
programmed  to  process  phantom  image  of  size  9  pels  by  9  pels  because  of  the 
limited  computer  main  memory  available  at  the  time.  With  our  recently 
enhanced  computer  facilities,  these  techniques  can  be  extended  to  a  larger  size 
of  sensor  array  that  would  occur  in  the  DSN  scenario. 

Relations  with  other  DSN  Contractors 

In  order  to  facilitate  experimentation  with  ladder  form  estimation  algo¬ 
rithms  on  actual  acoustical  sensor  data  and  to  take  advantage  of  existing  DSN 
software,  we  are  establishing  the  Lincoln  Laboratories  DSN  signal  processing 
software  at  Stanford.  Their  software  package  processes  acoustical  data 
recorded  on  digital  tape  from  their  microphone  arrays  during  aircraft  flyby 
experiments.  The  processing  of  data  from  a  single  sensor  array  produces  a 
power  versus  azimuth  plot  at  various  temporal  frequencies  and  elevations  of 
interest.  Appropriate  Interpretation  of  this  data  yields  bearing  information  for 
sound  sources  which  can  be  amalgamated  into  tracks  of  possible  target  motion. 
The  Lincoln  Lab  package  produces  only  bearing  information.  Since  the  max¬ 
imum  separation  between  sensors  in  the  current  array  is  three  meters,  it  is 
doubtful  that  high  resolution  range  information  for  distant  sources  could  also  be 
obtained  from  the  current  single  array  data  acquisition  system.  Using  the  data 
tape  reading  and  file  I/O  structure  of  this  package,  we  are  extending  the  power 
versus  azimuth  calculations  to  include  ladder  form  spectral  estimation  tech¬ 
niques. 

Our  efforts  to  establish  a  working  system  here,  first  on  a  VAX  and  later  on  a 
11/70  have  been  frustrated  by  many  difficulties  which  further  point  to  the 
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dty  o  establishing  a  common  signal  processing  workbench.  The  software 
[e  written  in  ”C”  for  a  11/70  running  UNIX  version  7  is  quite  large  and 
es  several  software  libraries  written  by  several  groups  at  Lincoln  Labs, 
torable  time  was  spent  in  locating  all  the  necessary  source  flies  and 
tly  transferring  them  over  the  ARPANET  initially  to  the  Mt.  St.  Helens  VAX 
iter.  The  location  of  the  various  source  components  is  listed  below.  The 
is  for  the  Tektronix  plotting  routines  came  from  the  Lincoln  Laboratories 
d  Seismology  group. 


ource 

eclarations  and  shell  flies 
irary  -list  management  system 
jrary  -parameter  data  base  subroutine 
•ary  -plotting  routine  sources 
from  11-asg 


/uO/green/dsnt 
/usr/bin 
/usr/src/lms 
/usr/src/pdbsubs 
/rml  3/crames/graph 
/gpactek 


be  attempt  to  compile  and  run  the  package  on  the  VAX  revealed  a  number 
tability  problems,  including  illegal  operations  on  pointers  and  incompati- 
e  of  data  types.  Manipulation  of  the  data  structures  and  file  header  infor- 
l  and  definition  of  file  I/O  variables  were  based  on  an  assumed  (11/70 
t)  2  byte  integer  word  length.  Since  the  VAX  word  length  is  4  bytes,  this 
loes  not  work  on  the  VAX  This  use  of  machine  dependent  code  appeared 
itly  so  that  conversions  were  substantially  more  involved  than  redeclaring 
le  types.  The  large  size  and  monolithic  construction  of  the  package  made 
nsequences  of  these  problems  severe.  Testing  for  program  bugs  was  com- 
sd  by  the  many  layers  of  subroutine  calls,  the  probability  of  multiple  bugs, 
ie  sheer  volume  of  code  to  search  Since  changes  to  one  procedure  tended 
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causes  a  need  for  modifications  in  other  procedures,  fixing  bugs  was  a  tedious 
>eration.  Even  when  changes  were  confined  to  one  procedure,  the  modification 
-ocess  was  slow  because  of  the  need  to  reload  the  entire  package.  AfUir  mak- 
g  a  list  of  the  problems  and  discussing  them  with  the  group  at  Lincoln  Labs,  we 
tached  an  agreement  that  they  would  rewrite  the  package  so  that  it  would  be 
unpiler  warning  free,  not  have  the  word  length  problems,  and  be  compatible 
ith  the  VAX  computer. 

The  recent  acquisition  of  a  11/70  computer  by  ISL  provided  the  opportunity 
j  attempt  to  bring  up  the  original  Lincoln  Laboratory  package  on  the  computer 
«■  which  it  was  written.  However  our  system  was  running  UNIX  version  6  rather 
lan  the  version  7  of  Lincoln  Labs.  We  experienced  problems  with  several 
lodules  that  were  in  the  version  7  operating  system  that  did  not  exist  in  version 
or  would  not  run  properly.  These  include  the  matrix  space  allocation  routine 
nd  the  ordered  library  building  routine.  The  large  number  of  source  files 
aused  problems  in  compiling  and  loading  the  program.  The  list  of  file  names 
xceeded  the  maximum  number  of  parameters  for  a  load  command,  so  libraries 
ad  to  be  constructed.  Due  to  the  large  program  size  and  the  restricted 
memory  of  the  DEC-1 1/70,  the  program  and  data  areas  have  to  be  separated. 
re  currently  have  the  azimuth  analysis  portion  of  the  Lincoln  Labs  software  Tun¬ 
ing  but  not  the  plotting  software.  Using  a  test  data  tape  of  artificial  data,  we 
tdnk  that  our  version  of  the  software  is  functioning  identically  to  the  Lincoln 
abs  software.  We  are  awaiting  a  data  tape  of  aircraft  flyby  experiments  in 
rder  to  apply  our  ladder  form  algorithms  on  real  sensor  data. 

Several  recommendations  can  be  made,  based  on  the  experience  of  tran- 
porting  this  package: 

Portability  should  be  a  high  priority  in  selection  of  programming  style.  This 

will  facilitate  both  transfer  of  software  between  groups  and  transfer  of 
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to  causes  a  need  for  modifications  in  other  procedures,  fixing  bugs  was  a  tedious 
operation.  Even  when  changes  were  confined  to  one  procedure,  the  modification 
process  was  slow  because  of  the  need  to  reload  the  entire  package.  After  mak¬ 
ing  a  list  of  the  problems  and  discussing  them  with  the  group  at  Lincoln  Labs,  we 
reached  an  agreement  that  they  would  rewrite  the  package  so  that  it  would  be 
compiler  warning  free,  not  have  the  word  length  problems,  and  be  compatible 
with  the  VAX  computer. 

The  recent  acquisition  of  a  11/70  computer  by  1SL  provided  the  opportunity 
to  attempt  to  bring  up  the  original  Lincoln  Laboratory  package  on  the  computer 
for  which  it  was  written.  However  our  system  was  running  UNIX  version  6  rather 
than  the  version  7  of  Lincoln  Labs.  We  experienced  problems  with  several 
modules  that  were  in  the  version  7  operating  system  that  did  not  exist  in  version 
6  or  would  not  run  properly.  These  include  the  matrix  space  allocation  routine 
and  the  ordered  library  building  routine.  The  large  number  of  source  files 
caused  problems  in  compiling  and  loading  the  program.  The  list  of  file  names 
exceeded  the  maximum  number  of  parameters  for  a  load  command,  so  libraries 
had  to  be  constructed.  Due  to  the  large  program  size  and  the  restricted 
memory  of  the  DEC-1 1/70,  the  program  and  data  areas  have  to  be  separated. 
We  currently  have  the  azimuth  analysis  portion  of  the  Lincoln  Labs  software  run¬ 
ning  but  not  the  plotting  software.  Using  a  test  data  tape  of  artificial  data,  we 
think  that  our  version  of  the  software  is  functioning  identically  to  the  Lincoln 
Labs  software.  We  are  awaiting  a  data  tape  of  aircraft  flyby  experiments  in 
order  to  apply  our  ladder  form  algorithms  on  real  sensor  data. 

Several  recommendations  can  be  made,  based  on  the  experience  of  tran¬ 
sporting  this  package: 

•  Portability  should  be  a  high  priority  in  selection  of  programming  style.  This 

will  facilitate  both  transfer  of  software  between  groups  and  transfer  of 
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software  to  new  machines  in  the  same  group. 

Installing  and  maintaining  software  is  easier  when  programs  are  small  and 
simple  in  structure. 

Complex  signal  processing  functions  should  be  achieved  through  linking  of 
simple  modules,  rather  than  by  building  large,  complex  programs. 

Last  but  not  least,  a  real-time  version  of  UNIX  should  be  developed. 
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6.  CONCLUSIONS 

We  have  made  considerable  progress  towards  the  development  of  new  algo¬ 
rithms  for  the  various  signal  processing  functions  acquired  in  a  DSN  system. 
Our  theoretical  research  on  the  source  location  problem  and  on  estimation 
techniques  has  reached  a  fairly  mature  stage.  We  now  need  to  turn  these  results 
into  software  modules  that  will  be  usable  by  Lincoln  Laboratory  and  other 
members  of  the  DSN  community.  The  signal  processing  workbench  will  enable 
us  to  test  the  algorithms  we  have  already  developed,  as  well  as  to  refine  them 
and  discover  new  ones.  We  need  to  evaluate  the  performance  of  these  algo¬ 
rithms  using  both  synthetic  and  real  data. 

Our  goal  is  to  gradually  replace  the  signal  processing  modules  currently 
used  by  Lincoln  Laboratory  with  our  own  improved  modules.  This  will  be  done  in 
stages.  First,  some  of  the  basic  components  will  be  replaced,  but  the  overall 
processing  structure  will  remain  unchanged.  For  example  a  high  resolution 
spectral  estimator  based  on  the  FFT  is  now  used  as  the  front  end  for  source 
location  detection.  We  can  replace  this  estimator  by  one  based  on  ladder  forms, 
which  is  expected  to  have  superior  performance.  In  a  second  stage  we  may  want 
to  make  more  comprehensive  changes,  which  may  require  structural 
modifications.  For  example,  a  combination  of  a  multichannel  ladder  form  with 
an  event  detection  algorithm  may  replace  the  current  detection  scheme.  A 
combination  of  an  ARMA  modeling  technique  for  TDOA  estimation  and  the  new 
algorithms  for  source  location  described  in  Section  D  may  replace  the  current 
localization  technique.  See  also  [DSN]  for  an  overview  of  the  signal  processing 
architecture  associated  with  these  new  algorithms.  This  process  of  replacing 
processing  modules  in  a  prototype  DSN  system  will  provide  a  conclusive  test  of 
our  algorithms,  and  lead  the  way  to  their  widespread  use  in  multi-source 
tracking /classification  applications.  A  test  of  this  type  is  an  essential  step  in 


the  development  of  any  new  technique. 

The  signal  processing  workbench  will  play  a  central  role  in  the  develop¬ 
ments  described  above.  It  will  enable  us  not  only  to  test  and  exercise  our  algo¬ 
rithms,  but  will  be  used  as  a  valuable  research  tool.  Some  of  the  algorithmic 
work  is  heavily  dependent  on  simulations  because  it  is  not  analytical.  Also, 
other  sophisticated  programming  languages  (e.g.,  MACSYMA)  will  be  needed  to 
handle  the  complexity  of  more  sophisticated  algorithms.  So  far  we  treated 
mainly  one  dimensional  problems.  The  development  of  two  and  three  dimen¬ 
sional  estimation  algorithms  is  expected  to  involve  complex  equation  manipula¬ 
tions  that  are  too  hard  to  perform  manually.  We  feel  that  our  progress  in  this 
area  will  be  limited  without  computer  aided  equation  handling. 

In  the  next  phase  of  our  project,  we  also  plan  to  emphasize  research  into 
some  of  the  fundamental  mathematical  issues  that  need  to  be  resolved  to 
achieve  advanced  DSN  system  design.  Much  of  the  earlier  work  on  multi-source, 
multi-sensor  signal  processing  involved  relatively  simple  modeling  (single  chan¬ 
nel,  stationary  processes,  etc.).  There  is  a  clear  need  to  develop  more  sophisti¬ 
cated  mathematical  models  and  the  analytical  tools  related  to  them.  Some 
examples  of  open  research  issues  are  briefly  mentioned  below: 

(i)  Performance  bounds  on  tracking  classification  algorithms,  given  multi¬ 
sensor  multi-source  information,  assuming  non-stationary  and  non-gaussian 
processes.  Such  bounds  are  essential  for  the  evaluation  of  the  performance 
of  any  DSN  system,  i.e„  how  far  away  from  "optimal"  is  it  ? 

(ii)  Additional  insights  are  needed  into  the  question  of  distributed  processing. 
We  need  to  find  ways  of  characterizing  the  minimal  information  that  needs 
to  be  exchanged  in  a  distributed  system  in  order  to  achieve  a  given  perfor¬ 
mance  level.  How  much  of  the  information  can  be  used  only  at  the  local 
level,  and  how  much  needs  to  be  communicated  to  other  nodes. 
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(iii)  Development  of  efficient  estimation  algorithms  for  two  and  three  dimen¬ 
sional  systems.  We  need  to  find  a  good  spectral  representation  for  multi¬ 
dimensional  processes.  It  seems  clear  that  there  are  differences  at  a  very 
fundamental  level  between  the  structure  of  multi-dimensional  and  one 
dimensional  problems.  As  mentioned  before,  the  complexity  of  the  equa¬ 
tions  arising  in  such  systems  makes  it  difficult  to  search  for  structural  pat¬ 
terns. 

(iv)  A  different  approach  to  distributed  computations,  that  are  not  explored  so 
far.  is  to  use  probabilistic  algorithms  of  the  type  developed  by  Rabin. 
Independent  solutions  by  multiple  nodes  followed  by  majority  voting  can 
lead  to  a  natural  way  of  distributing  computations  in  a  complex  algorithm. 

(v)  A  large  class  of  DSN  problems  combines  both  processing  and  communica¬ 
tion  of  data.  The  communication  aspects  of  these  problems  involve  the 
theory  of  multiterminal  communication.  The  main  results  available  so  far 
show  bounds  on  the  channel  capacity  (the  maximum  rate  of  information 
transmission)  for  broadcast  or  multi-access-type  channels.  Most  of  these 
results  are  of  the  existence  type,  l.e.,  do  not  provide  practical  schemes  that 
could  be  implemented.  See  [VM]  [BER]  [AEG]. 

In  summary,  we  plan  to  continue  our  investigation  on  two  levels,  basic 
research  into  fundamental  mathematical  issues  of  the  DSN,  and  the 
development/implementation/testing  of  signal  processing  modules  for  the  pro¬ 
totype  DSN  system.  This  approach  combines  long  term  goals  such  as  optimal 
DSN  design  methodology  with  short  term  performance  gains. 


•  "•*  ,*  "  .*  *  '  •  ' .  *  “  m  *  »  *  •  *  •  *  *  *  •  *  •  *  •  *  •  *  •  *  •  *  •  **  *  *'  ■  '  *  *  »'•  **  i**  »*•  .**  .**  •"'*  -  ’>  *  .*•>  *,*  i 

4  *  •  *  .  *  .  *  •  *  •  .  *  ■  'k  .*.'«***»•»*  I*  J  ■  >>  k  *  |  *  ,  *  .  »  .V  •  *  1 


-96- 


7.  DSN  PUBLICATIONS 

This  section  summarizes  the  publications  and  presentations  resulting  from 
our  work  on  the  DSN  project.  With  few  exceptions,  these  publications  are 
described  in  Sections  2,  3  and  4. 


SUBMITTED: 

B.  Porat,  M.  Morf  and  D.  Morgan.  "On  the  Relationship  Among  Several  Normalized 
Square-Root  Ladder  Algorithms,"  CISS,  submitted  January  1981. 

A.  Nehorai  and  M.  Morf,  "Enhancement  of  Sinusoids  in  Colored  Noise  and  Whiten¬ 
ing  Performance  of  Least-Squares  Predictors,"  IEEE  Trans.  ASSP,  submitted 
January  1981. 

B.  Porat,  B.  Friedlander  and  M.  Morf,  "Square  Root  Covariance  Ladder  Algo¬ 
rithms."  IEEE  Trans  on  Automatic  Control,  submitted  December  1980. 

B.  Egardt  and  M.  Morf.  "Asymptotic  Analysis  of  A  Ladder  Algorithm  for  ARMA 
Models,"  IEEE  Trans,  on  Automatic  Control,  submitted  November  1980. 

M.  Morf,  "Displacement  Ranks  and  Doubling  Algorithms  for  Matrix  Inversion," 
Bulletin  of  the  American  Mathematical  Society,  submitted  1979. 

M.  Morf.  "Doubling  Algorithms  for  Toeplitz  and  Related  Equations,"  IEEE  Trans. 
Inf.  Theory ,  to  appear  1981. 

R.  Longchamp  and  M.  Morf,  "Distributed  Gaussian  Second-Order  filters,"  IEEE 
TYans.  on  Automatic  Control,  submitted  December  1979. 

M.  Morf.  "Extended  System  and  Transfer  Function  Matrices  and  System 
Equivalence,"  Int.  J.  of  Control,  submitted  November  1979. 

M.  Morf,  ’Tast  Cholesky  Algorithm  and  Adaptive  Feedback  Ladder  Forms,"  IEEE 
Proc.' 80th  CDC.  Dec.  1981,  submitted  March  1981. 

P.R.R.L.  Nunes  and  M.  Morf,  "Image  Reconstruction  Techniques  and  Target  Detec¬ 
tion.”  IEEE  Pro c.  80th  CDC,  Dec.  1981,  submitted  March  1981. 

D.T.L  Lee,  J.-M.  Delosme  and  M.  Morf,  "State-Space  Structures  of  Generalized 
Ladder  Canonical  Forms,"  IEEE  Proc.  20th  CDC,  Dec.  1981,  submitted  March 
1981. 

J.-M.  Delosme  and  M.  Morf,  "Normalized  Doubling  Algorithms  for  Finite  Rank 
Processes,"  IEEE  Proc.  80th  CDC,  Dec.  1981,  submitted  March  1981. 

H.  M.  Ahmed  and  M.  Morf,  "On  the  Relationship  Among  Certain  Numerical  Algo¬ 
rithm  and  Bilinear  Bang-Bang  Control,"  IEEE  Proc.  80th  CDC,  Dec.  1981,  submit¬ 
ted  March  1981. 
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M.  Morf,  C.H.  Muravchik,  D.T.  Lee  and  J.  M.  Delosme,  "Hilbert  Space  Array 
Methods  for  Finite  Rank  Process  Modeling  and  Ladder  Form  Realizatins."  IEEE 
Proc.  20th  CDC  and  IEEE  Trans.  Automatic  Control,  Dec.  1981,  submitted  March 
1981. 

M.  Morf  and  A.  Nehorai,  "A  Nev  Algorithm  for  Recursive  Estimation  of  Right 
Matrix  Fraction  Description  Type  Multivariable  Systems."  IEEE  Pro c.  20th  CDC, 
Dec.  1981,  submitted  March  1981. 

Mohamed  T.  Hadidi  and  M.  Morf,  "Efficient  Solution  of  the  Inverse  Levinson  Prob¬ 
lem  in  the  Multichannel  Case  Via  the  Lyapunov  Equation,"  IEEE  Proc.  20th  CDC, 
Dec.  1981,  submitted  March  1981. 


ACCEPTED: 

H.M.  Ahmed,  P.H.  Ang  and  M.  Morf,  "A  VLSI  Speech  Analysis  Chip  Set  Utilizing  Co¬ 
ordinate  Rotation  Arithmetic,"  ICCAS,  Chicago,  IL.  April  1981.  (invited) 

B.  Porat,  B.  Friedlander  and  M.  Morf,  "Square  Root  Covariance  Ladder  Algo¬ 
rithms."  Proc.  1981  ICASSP,  Atlanta,  GA,  March  30,  1981,  pp.  877-880. 

B.  FYiedlander  end  M.  Morf,  "Efficient  Algorithms  for  Adaptive  Linear-Phase 
Filtering,"  Proc.  1981  ICASSP,  Atlanta.  GA.  March  30.  1981,  pp.  247-250. 

A.  Nehorai  and  M.  Morf.  "Enhancement  of  Sinusoids  in  Colored  Noise  and  Whiten¬ 
ing  Performance  of  Least-Squares  Predictors,"  Proc.  1981  ICASSP,  Atlanta,  GA, 
March  30.  1981.  pp.  275-278. 

M.  Morf,  C.  Muravchik  and  D.T.  Lee.  "Hilbert  Space  Array  Methods  for  Alpha- 
Stationary  Process  Estimation  and  Ladder  Realizations  for  Speech  and  Adaptive 
Signal  Processing,"  Proc.  1981  ICASSP,  Atlanta,  GA.  March  30,  1981,  pp.  B56-B59. 

H.M.  Ahmed.  M.  Morf,  D.T.  Lee  and  P.  Ang.  "A  VLSI  Speech  Analysis  Chip  Set  Based 
on  Square-Root  Normalized  Ladder  Forms,”  Proc.  1981  ICASSP,  Atlanta.  GA, 
March  30.  1981,  pp.  648-653. 

D.T.  Lee  and  M.  Morf,  "Alpha-Stationary  Distortion  Measures  via  Ladder  Forms," 
IEEEInt.  Sym.  on  Inf.  Theory,  Feb.  9-12,  1981. 

M.  Morf  and  J.M.  Delosme,  "Pseudo  Innovations  Representations  for  Alpha- 
Stationary  Processes,"  IEEE  Int.  Sym.  on  Inf.  Theory,  Feb.  9-12,  1981. 

D. T.L  Lee.  M.  Morf  and  B.  FYiedlander,  "Recursive  Square-Root  Ladder  Estima¬ 
tion  Algorithms,"  special  joint  issue  of  IEEE  Trans.  C&S  and  Trans.  ASSP,  to 
appear  in  June  1981. 

S.  Wood  and  M.  Morf,  "A  Fast  Implementation  of  a  Minimum  Variance  Estimator 
for  Computerized  Tomography  Image  Reconstruction"  IEEE  Trans,  on  Bio- 
Medical  Engineering,  Vol.  BME-28,  No.  2,  pp.  56-68,  February  1981. 

E.  Verriest,  B.  Friedlander  and  M.  Morf,  "Distributed  Processing  in  Estimation 
and  Detection."  18th  IEEE  Conference  on  D&C,  Dec.  1979. 


J.M.  Delosme  and  M.  Morf,  "Mixed  and  Minimal  Representations  for  Toeplitz  and 
Related  Systems,"  Proc.  14th  Asilomar  Con/.  on  CSC,  Nov.  17-19,  1980. 

J.-M.  Delosme,  Y.  Genin,  M.  Morf  and  P.  Van  Dooren,  "Sigma-Contractive  Embed¬ 
dings  and  Interpretation  of  Some  Algorithms  for  Recursive  Estimation,"  14th  Asi¬ 
lomar  Can/,  on  Circuits,  Systems  and  Computers,  Nov.  17-19,  1980. 

S.L.  Wood  and  M.  Morf,  "A  Fast  Implementation  of  the  Minimum  Variance  Estima¬ 
tor  for  CT  Image  Reconstruction  Application,"  Proc.  14th  Asilomar  Con/,  on  CSC, 
Nov.  17-19.  1980. 

B.  Friedlander  and  M.  Morf,  "Efficient  Inversion  Formulas  for  Sums  of  Products  of 
Toeplitz  and  Hankel  Matrices,"  18th  Annual  Allerton  Conference.  October  8-10, 
1980. 

H.M.  Ahmed,  J.-M.  Delosme  and  M.  Morf,  "Highly  Concurrent  Computing  Struc¬ 
tures  for  Digital  Signal  Processing  and  Matrix  Arithmetic,"  IEEE  Computer, 
March  1981.  (invited) 


PUBLISHED 

D.  Lee.  B.  Friedlander  and  M.  Morf.  "Recursive  Ladder  Algorithms  for  ARMA 
Modeling."  19th  Proc.  IEEE  CDC,  December  10,  1980,  pp  1225-1231. 

D.  Lee  and  M.  Morf.  "State-Space  Structures  of  Ladder  Canonical  Forms,"  19th 
Proc.  IEEE  CDC,  December  10.  1980,  pp  1221-1224. 

M.  Morf  and  J.M.  Delosme,  "A  Tree  Classification  of  Algorithms  for  Toeplitz  and 
Related  Equations  Including  Generalized  Levinson  and  Doubling  Type  Algo¬ 
rithms."  19th  Proc.  IEEE  CDC,  December  10.  1980,  pp  42-46. 

J.  Turner,  "Use  of  the  Digital  Lattice  Structure  in  Estimation  and  Filtering." 
Proc.  European  Signal  Processing  Con/.,  Switzerland,  Sept.  1980,  pp.  33-42. 

R.  Longchamp.  "Stable  Feedback  Control  of  Bilinear  Systems,”  IEEE  Trans,  o/ 
Automatic  Control,  Vol.  AC-25,  No.  2,  pp.  302-306,  April  1980. 

R.  Longchamp,  "Controller  Design  for  Bilinear  Systems,"  IEEE  Automatic  Con¬ 
trol,  in  Vol.  25.  No.  3,  June  1980,  pp.  547-548. 

M.  Morf.  "Doubling  Algorithms  for  Toeplitz  and  Related  Equations"  Proc.  1980 
ICASSP,  Denver.  CO,  April  9-11,  1980,  pp  954-959. 

D.T.L  Lee  and  M.  Morf,  "A  Novel  Innovations  Based  Time-Domain  Pitch  Detection" 
Proc.  1980  ICASSP.  Denver.  CO,  April  9-11,  1980,  pp.  40-44. 

J.M.  Delosme.  M.  Morf  and  B.  Friedlander,  "Source  Location  from  Time  Difference 
on  Arrival:  Identiflability  and  Estimation"  Proc.  1980  ICASSP,  Denver,  CO,  April 
9-11,  1980,  pp.  818-824. 

M.  Morf  and  D.T.L  Lee,  "Recursive  Square-Root  Ladder  Estimation  Algorithms," 
Pro c.  1980  ICASSP,  Denver,  CO,  April  9-11,  1980,  pp.  1005-1017. 
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M.  Morf  and  D.  Lee,  "Recursive  Least  Squares  Ladder  Forms  for  Fast  Parameter 
Tracking,"  Proc.  IEEE  Con f.  D&C,  San  Diego,  CA,  Jan.  10-12,  1979,  pp.  1362-1367. 

M.  Morf  and  J.M.  Delosme,  "Doubling  Algorithms  for  Block  Toeplitz  and  Related 
Equations  with  Applications  to  2-D  Problems,"  Presented  Workshop  on  2-D 
Dimensional  Digitial  Signal  Processing  at  Lawrence  Hall  of  Science,  Oct.  1979. 

R.  Longchamp  and  M.  Morf,  "Estimators  for  Quadratic  Dynamic  Systems," 
Presented  13th  Annual  Asilomar  Conf.,  Nov.  1979,  pp.  463-468. 

M.  Morf  and  D.  Lee.  "Recursive  Spectral  Estimation  of  Alpha-  Stationary 
Processes,”  Proc.  Rome  Air  Development  Center  Spectrum  Estimation 
Workshop,  pp.  97-108,  1978. 

M.  Morf,  B.  Friedlander  and  J.  Newkirk,  "Tutorial  Survey  of  Algorithms  for  Locat¬ 
ing  and  Identifying  Spatially  Distributed  Sources  and  Receivers,"  ARPA  Distri¬ 
buted  Sensor  Network  Symposium,  pp.  94-104,  1978. 

M.  Morf,  D.  Lee,  J.  Nickolls  and  A.  Vieira,  "A  Classification  of  Algorithms  for  ARMA 
Models  and  Ladder  Realizations."  Proc.  IEEE  ICASSP.  pp.  13-19,  May  1977. 
Reprinted  in  IEEE  Press  Series.  Modem  Spectrum  Analysis  edited  by  D.G.  Child¬ 
ers,  pp.  262-268,  1978. 


THESES: 

Levy.  B..  "Algebraic  Approach  to  Multi-Dimensional  Systems."  June  1981. 

Lee.  D.T.L.,  "Ladder  Form  Realizations  of  Fast  Algorithms  in  Estimation,"  August 
1980. 

Verriest.  E.,  "Structure  of  Time  Varying  Systems,"  August  19B0. 

Fortes.  J.M.P.,  "An  Estimation  Approach  to  3-D  Reconstruction  Problems  Includ¬ 
ing  Counting  Statistics  With  Applications  to  Medical  Imaging,"  June  I960. 

Nunes,  P.R.,  "Estimation  Algorithms  for  Medical  Imaging,"  June  1980. 

J.  Newkirk.  "Computational  Issues  In  Least-Squares  Estimation  and  Control." 
June  1979. 


REPORTS: 

M.  Morf,  B.  Friedlander  and  J.  Newkirk,  'Investigation  of  New  Algorithms  for 
Locating  and  Identifying  Spatially  Distributed  Sources  and  Receivers,"  Annual 
Technical  Summary  Report  to  DARPA,  SEL  Report  No.  M355-1,  March  31,  1979. 

Debate ur.  S..  "An  Analysis  of  Spatial-Temporal  Filtering  in  Remote  Probing  of 
Atmospheric  Turbulence  and  Wind  Speed,"  Technical  Report  4510-1,  Stanford 
Electronics  Lab,  Stanford  University,  December  1980. 
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